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ABSTRACT 


In  recent  years  the  newly  developed  Singularity  Expansion  Method 
(SEM)  has  been  applied  numerically  to  many  transient  wire  antenna  and 
scattering  problems.  The  present  work  applies  it  to  a two-dimensional 
problem  - that  of  the  electromagnetic  penetration  through  a rectan- 
gular aperture. 

The  aperture  problem  is  approached  through  the  electromagnetic 
dual  structure  - the  rectangular  plate.  The  SEM  quantities  for  it  are 
determined  numerically  by  way  of  the  Integrated  form  of  the  electric 
field  Integral  equations.  Numerical  efficiency  considerations  are 
important  in  the  solution  of  a two-dimensional  problem.  The  important 
details  relating  to  computational  efficiency  that  were  applied  in  the 
present  work  are  described. 

The  pole  trajectories  for  the  plate  are  given  for  a range  of  aspect 
ratios,  along  with  natural  modes  at  specimen  poles.  Physical  and  nu- 
merical interpretations  for  these  data  are  stated  where  possible. 

Theoretical  aspects  of  SQ1  relating  to  coupling  coefficient  forms 
and  to  the  role  of  complex  domain  entire  functions  in  the  time  domain 
SEM  representation  are  explored.  Correct  coupling  coefficient  forms 
are  found  to  depend  on  the  asymptotic  behavior  at  infinity  of  the  fre- 
quency domain  Green's  function.  This  asymptotic  behavior  is  unrecoverable 
in  numerical  solutions,  however.  Tests  on  the  time  domain  expansion  of 


the  Green's  function  are  proposed  herein  to  determine  the  validity  of  a 
Par^cu^-ar  coupling  coefficient  form  for  a particular  problem  when  the 
asymptotic  behavior  is  unavailable.  It  is  observed  that  Baum's  class  1 
coupling  coefficient  can  be  made  valid  in  problems  where  it  otherwise 
would  not  apply  at  the  expense  of  including  an  entire  function  contribu- 
tion in  the  time  domain  expansion.  This  entire  function  contribution  is, 
in  general,  significant. 
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1.  INTRODUCTION 

1 . 1 Description  of  the  Work 

This  document  presents  the  results  of  an  investigation  of  repre- 
senting the  transient  electromagnetic  coupling  through  a rectangular 
aperture  in  a conducting  screen  by  means  of  the  Singularity  Expansion 
Method  [1].  The  problem  is  formulated  by  means  of  its  electromagnetic 
dual  problem,  the  scattering  from  a rectangular  plate.  Consequently, 
the  results  are  equally  applicable  to  both  members  of  the  dual  pair. 

The  Singularity  Expansion  Method  (SEM)  was  introduced  by  Baum  in 
1971  [2].  It  provides  a compact  means  of  representing  broadband  transient 
electromagnetic  phenomena  on  resonant  bodies  in  terms  of  the  complex 
natural  resonances  of  the  body  in  the  complex  Laplace  transform  plane. 

For  the  problem  at  hand,  these  resonances  are  located  numerically  by 
means  of  a method  of  moments  solution  to  an  integral  equation  formula- 
tion. The  integral  equations  used  to  model  the  plate  structure  are  the 
integrated  forms  of  the  coupled  electric  field  integral  equations  for  the 
plate  developed  by  Rahmat-Samii  and  Mittra  [3].  Their  original  formulation 
holds  for  the  real  frequency  domain,  but  extends  directly  to  complex 
frequencies  for  the  purposes  herein. 

1 . 2 The  Singularity  Expansion  Method 

Subsequent  to  its  introduction  by  Baum  in  1971,  SEM  has  been  applied 
to  many  problems  - both  numerically  and  analytically.  In  fact,  in 
his  introductory  paper  [2],  Baum  constructed  the  formal  SEM  solution  for 
a perfectly  conducting  spherical  scatterer.  Marin  conducted  an  analytical 
solution  for  a prolate  spheroidal  scatterer  [4].  Umashankar  has  applied 
SEM  analytically  to  the  circular-loop  wire  antenna  [5]. 
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Marin  and  Latham  [6]  made  an  important  theoretical  contribution  to 
SEM  early  in  its  development  . They  demonstrated  that  in  the  complex 
frequency  plane  the  only  singularities  in  fields  for  finite  perfectly 
conducting  bodies  are  poles.  This  fact  precludes  the  need  to  consider 
branch-cut  integrals  and  allows  broader  applications  of  analytic  con- 
tinuation concepts  in  the  frequency  plane. 

Tesche  performed  the  first  numerical  SEM  construction  [7].  He 
considered  the  problem  of  a thin  cylindrical  antenna.  Subsequently, 
various  workers  have  applied  the  method  numerically  to  the  L-wire  [5], 

16],  [9],  and  several  coupled  cylinder  problems  [ 10 ] — [ 12 ] . Tesche  ex- 
plored SE>1  solutions  to  resistively  loaded  thin  wires  [13]. 

In  previous  work,  surface  problems  have  been  treated  by  analytical 
methods, and  numerical  methods  have  been  applied  only  to  wire  structures. 

One  contribution  of  the  present  work  is  to  indicate  the  viability  of 
the  numerical  approach  to  SEM  on  a surface  problem. 

1 . 3 Scope  of  the  Present  Work 

There  are  three  chapters  of  original  substance  in  this  document. 

Chapter  2 delves  into  some  theoretical  considerations  that  are  important 
to  all  SEM  problems.  Chapters  3 and  4 deal  with  the  specific  problem 
of  the  rectangular  aperture. 

Chapter  2 defines  the  SEM  formalism  and  presents  some  new  insights  into 
the  inversion  procedure.  In  particular,  it  relates  the  asymptotic  behavior  of 
the  frequency  domain  SEM  representation  to  the  definition  of  the  so-called 
coupling  coefficient.  Because  all  of  the  terms  which  enter  into  the 
asymptotic  form  are  not  recoverable  in  numerically  derived  SEM  solutions,  we 
must  turn  to  indirect  means  of  defining  the  role  of  the  asymptotic  behavior. 
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It  is  shown  in  Chapter  2 that  the  asymptotic  behavior  of  the  frequency 
domain  representation  dictates  the  support  of  the  time  domain  Green's 
function  for  a given  problem. 

We  introduce  two  means  of  overcoming  the  difficulties  presented  by 
the  unavailable  asymptotic  information.  Cne  of  the  means  is  a postulated 
asymptotic  behavior  which  is  conservative  and  physically  consistent.  This 
postulated  form  results  in  a new  "class  X"  coupling  coefficient  which  may 
be  used  directly  or  as  a basis  of  comparison  for  other  coupling  coefficient 
forms.  The  second  means  is  a test  on  the  time  domain  expansion  to  deter- 
mine for  a specific  problem  whether  or  not  any  particular  coupling  co- 
efficient form  is  consistent  with  the  physical  concept  of  causality.  The 
class  X coupling  coefficient,  the  effect,  imposes  causality  directly. 

Chapter  3 presents  the  complex  frequency  generalization  of  the 
thin  plate  integral  equations  and  describes  the  numerical  procedures 
used  in  determining  the  SEM  quantities  from  this  formulation.  Because 
the  present  problem  is  a surface  problem,  we  are  forced  to  be  particularly 
conscious  of  numerical  efficiency.  Some  algorithmic  considerations 
that  are  important  to  numerical  efficiency  are  discussed  in  Chapter  3. 

The  argument  number  procedure  proposed  by  Baum  [14]  was  used  to  advantage 
in  isolating  the  SEM  poles  in  the  course  of  this  work.  The  details  of  its 
application  are  given  in  Section  3.5. 

Chapter  4 contains  numerically  obtained  SEM  data  for  the  rectangular 
plate.  Complete  pole  data  for  a portion  of  the  complex  plane  are  given 
for  rectangles  with  aspect  ratios  between  0.1  and  1.0.  Representative 
natural  mode  data  are  given  to  provide  the  reader  a means  of  physical 
Interpretation  of  the  resonances  of  the  structure.  Complete  mode  data 
and  normalization  constants  are  available  for  the  plate  in  computer  data 
base  form  [15]. 
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1.4  Notation 

As  in  many  other  fields  problem  areas,  mathematical  notation  for 
SEM  that  is  simultaneously  compact  and  explicit  is  desirable  but  a bit 
elusive.  The  notation  used  herein  is  built  on  that  of  Baum  [e.g.,  2]. 
Symbols  denoting  functions  are  chosen  for  their  connotative  value  where 
possible  and  their  attributes  are  indicated  with  bars,  tildes,  etc. 

In  particular,  a vector  quantity  is  indicated  with  a single  overbar, 
and  a dyadic  quantity  by  two  overbars.  Laplace  transform  pairs  are 
denoted  with  the  same  alphabetic  and  vector  attribute  symbols  in  both 
domains.  A tilde  (~)  oversymbol  identifies  the  function  of  the  pair 
which  has  as  its  domain  complex  frequency,  while  the  absence  of  the 
oversymbol  denotes  its  time  domain  counterpart. 

The  inner  product  notation  ^A(r),  B(r)^  means,  per  Baum,  to  take 
the  vector  inner  product  and  integrate  over  the  common  spatial  support. 

We  have  occasion  to  indicate  explicitly  time  domain  convolution  and 
complex  contour  integration  in  addition  to  the  spatial  inner  product. 

For  example, 

QD 

E(r,t)  = / (r(r,r’,t  - t’),  J(r’,t)^  dt ' 

—00 

indicates  that  the  time  domain  vector  quantity  E(r,t)  is  the  time  domain 

ao 

convolution  / dt’  of  the  vector  dot  product  of  the  dyad  F and  the 

—ao 

vector  J integrated  over  the  spatial  support  common  to  J and  F.  The 
transform  domain  counterpart  to  this  expression  is 

E(r,s)  - (f(r,r',s),  j(r’,s>) 

where  the  tildes  indicate  transform  domain  image  functions.  The  convolu- 
tion in  time  appears  as  a simple  multiplication  in  the  complex  frequency 


domain. 
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2.  THE  SEM  REPRESENTATION  AND  INVERSION  PROCEDURE 
2 . 1 Int  roduct ion 

In  this  chapter,  we  develop  the  Laplace  transform  solution  for  a 
general  linear  operator  equation  involving  time  and  spatial  variations. 

The  point  of  view  taken  here  departs  from  those  of  the  previous  workers 
addressing  SEM  solutions.  In  any  SEM  development,  the  primary  concern 
is  to  construct  the  solution  to  the  operator  equation  in  a Laplace  trans- 
form domain,  where  the  transform  is  taken  with  respect  to  the  time  vari- 
able in  the  equation.  Where  previous  work  has  viewed  this  transform 
solution  as  the  composition  of  a residue  series  and  an  entire  function 
in  the  transform  variable,  the  present  work  views  the  transform  domain 
solution  as  the  product  of  a driving  term  and  the  transform  of  a causal 
Green's  function  or  inverse  kernel  (see,  for  example,  Stakgold  [1]). 

The  Green's  function  may  be  formally  constructed  as  an  eigenfunction 
expansion  [2].  The  inverse  transform  to  the  time  domain  may  be  taken 
by  appealing  to  Jordan's  lemma  and  the  Cauchy  residue  theorem. 

Indeed,  the  last  step  introduces  one  of  the  most  significant  features 
in  the  SEM  approach.  The  appeal  to  Jordan's  lemma  [3]  depends  on  the 
asymptotic  form  of  the  transform  of  the  solution  for  large  s.  The  causal 
behavior  of  the  time  domain  solution  is  imposed  through  this  asymptotic 
form. 

The  asymptotic  behavior  of  the  Green's  function  factor  in  the  solu- 
tion can,  in  principle,  arise  either  from  the  infinite  pole  series  or 
from  the  entire  function  term  in  the  Mittag-Lef fler  representation  [A] 
of  the  Green's  function.  For  example,  the  Green's  function  for  the 
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wire  loop  shows  exponential  asymptotic  behavior  resulting  from  the  pole 
series  [5].  For  problems  which  are  solved  numerically,  only  a finite 
number  of  poles  in  this  representation  are  available  to  us.  Consequently, 
the  asymptotic  behavior  of  the  transform  solution  is  not  directly  recover- 
able. In  Section  2.5.2  we  develop  heuristically  the  dominant  asymptotic 
form  for  the  Green’s  function  of  an  integral  operator  with  a free-space 
wave  equation  related  kernel.  This  asymptotic  form  leads  to  a new  type 
of  coupling  coefficient  for  the  SEM  representation.  This  new  "class  X" 
coupling  coefficient  is  observed  to  impose  causality  directly  on  the 
solution.  We  observe,  however,  that  other  coupling  coefficient  forms 
that  do  not  impose  causality  directly  have  been  used  with  varying  success 
by  other  workers.  A procedure  determining  when  the  new  form  may  be  re- 
duced to  a simpler  coupling  coefficient  form  is  presented.  The  pro- 
cedure is  based  on  observation  of  the  time  domain  transform. 


2. 2 Transform  Solution  of  Time  Domain  Electromagnetic  Integral  Equations 
2.2.1  Laplace  transform  of  integral  equation 

We  state  a general  time  domain  integral  equation  for  perfectly  con- 
ducting finite  extent  bodies  in  free  space  as  follows. 

CO 

/ (f(r,r',t  - t'),  j(r',t’)^  dt ' = E(r,t)  , r e B, 

— no 


f(r,r\t) 

J(r,t) 

E(r,t) 


= (Tl 
= 0 \ 


= 0 


r,r'  e B and  t < 0 


(2.1) 
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The  brackets  denote  an  integration  of  the  vector  inner  product  of  the  two 
operands  over  the  support  of  the  second  operand.  The  P term  is  a dyadic 
kernel  to  the  integral  equation;  J is  a characteristic  electromagnetic 
quantity,  generally  a current  density;  and  E is  an  incident  wave  or 

CD 

driving  function  field.  The  J dt'  is  a time  domain  convolution.  The 

— QO 

body  is  denoted  by  the  region  B in  r space. 

We  proceed  to  construct  the  inversion  of  (2.1).  Provided  the 
Laplace  transforms  exist  for  both  the  left  and  right  sides  of  the  in- 
tegral equation,  we  may  write  it  in  the  transform  domain 

\r(r,r',s),  J(r',s)^  = E(r  ,s)  , r e B , Re  s > 0 . (2.2) 

The  tilde  (~)  overbar  is  used  to  denote  the  Laplace  transform  domain 
images  of  the  associated  functions  from  (2.1).  In  the  familiar  way,  the 
convolution  over  the  time  variable  goes  over  to  an  algebraic  product 
in  the  transform  domain. 

We  appeal  to  the  physical  concept  of  causality  and  observe  that 
each  of  the  three  functions  in  question  must  be  causal  in  nature  so  long 
as  we  require  that  E(r,t)  be  zero  for  finite  r and  t less  than  some 
finite  starting  time.  If  all  of  our  functions  are  causal,  we  may  with 
an  appropriate  shift  of  origin  convert  our  transform  to  a one-sided 
one.  Therefore,  it  is  sufficient  to  consider  all  transforms  as  one- 
sided ones,  and  we  define,  explicitly, 

00 

{ f (t ) } - f(s)  = / f(t)  e_st  dt  , Re  s > 0 (2.3) 

0 

under  the  condition  that  f(t)  be  bounded  for  t > 0. 
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2.2.2  An  inversion  of  the  transformed  Integral  equation 

We  may  write  a formal  solution  to  (2.2)  as 

J (r,s)  = \F  1(r,r',s),  E(r',s))  (2.4) 

where  the  inverse  kernel  T ^ (Green's  function)  is  the  solution  to 

( ^(r.r'  >s),  r(r',r0,s))  = I6(r  - rQ)  , (2.5) 

I is  an  identity  dyad  on  the  three  spatial  vectors  and  6(r  - r^)  is  a 
three-dimensional  Dirac  distribution.  The  equations  (2.4)  and  (2.5)  hold 
in  a distributional  sense.  In  fact,  Schwartz  points  out  that  F ^ cannot 
be  a function  [6].  In  addition,  Marin  and  Latham  have  shown  that  for 
perfectly  conducting  bodies  of  finite  extent,  T (r,r',s)  is  analytic 
in  s except  for  pole  singularities  [7].  This  latter  property  gives  us 
the  means  of  expressing  our  desired  result,  the  time  domain  solution  to 
(2.1)  in  terms  of  residues  at  the  poles  of  T . Since  the  body  is  pas- 
sive and  since  a pole  gives  rise  to  temporal  behavior  exp(s^t),  we 

conclude  that  all  of  the  poles  of  F ^ have  a negative  real  part  so  as 
to  produce  decaying  time  behavior.  Further,  since  (2.1)  has  a real 
solution,  it  is  evident  that  the  poles  must  occur  in  conjugate  pairs  in 
the  complex  s plane. 

For  many  problems  of  practical  interest,  we  can  recover  by  numerical 
or  analytical  methods  the  poles  which  dominate  the  time  response  for 
excitations  whose  spectra  are  bandlimited.  The  key  to  practical  appli- 
cation of  SEM  is  to  construct  the  time  domain  solution  to  (2.1)  in  such 
a way  that  it  is  accurately  represented  by  this  hopefully  modest  collection 


of  poles. 


2.2.3  The  Laplace  inversion  of  the  solution 


The  time  domain  solution  to  (2.1)  is  constructed  by  forming  the 
Laplace  inversion  of  (2. A)  as  follows. 

J(r,t)  - / J(r,s)  eSt  ds 

CB 

= 1(r»r',s),  E(r',s))  eSt  ds  . (2. 

CB 

The  integration  is  over  the  Bromwich  contour,  C , given  by 

B 

s = a + ju 

a > 0 , u e (-00,00) 

The  real  value  of  s for  this  contour  must  be  chosen  so  as  to  place  CD 
in  the  region  of  convergence  for  the  representation  of  J(r,s).  Conse- 
quently, we  must  have  a representation  for  J(r,s)  which  is  valid  in 
the  right-half  plane,  i.e.,  that  is  convergent  and  solves  (2.2)  for  s 
such  that  a > Oq,  or  a fixed  positive  number. 

Typically,  in  the  Laplace  transform  theory,  one  appeals  to  an  in- 
direct evaluation  of  the  integral  of  (2.6).  Namely,  C_  is  closed  with 
a semicircular  arc  of  infinite  radius  either  to  the  right  or  to  the  left 
as  indicated  in  Figure  2.1.  The  asymptotic  behavior  of  the  integrand 
generally  dictates  that  the  integrand  vanish  on  one  or  the  other  of  the 
semicircles  so  that  Jordan's  lemma  [3]  is  satisfied  and  the  integration 
over  either  the  closed  contour  or  Cg  + Cm  is  equivalent  to  in- 
tegration over  Cg  alone.  The  composite  contours  and 


-plane 
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The  Bromwich  contour  for  inversion  of  the  Laplace  transform  and  possible  closures  thereof. 
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are  closed  contours.  Hence,  when  Integration  along  C is  equivalent  to 

D 

integration  over  one  of  these  closed  contours,  we  may  evaluate  the  in- 
tegration by  means  of  Cauchy's  residue  theorem  so  that  information  at  the 
poles  of  J(r,s)  suffices  to  represent  J(r,t). 

Before  considering  the  asymptotic  form  of  J(r,s),  we  shall  develop 
the  singularity  expansion  in  a formal  manner.  The  eigenfunction  ex- 
pansion point  of  view  [2]  is  taken  in  this  development. 

The  formal  solution  to  (2.2)  may  be  constructed  by  means  of  an 
eigenfunction  expansion.  That  is,  we  define  the  left-  and  right-hand 
eigenvectors  L^Cr.s),  R^(r,s)  by  the  eigenvalue  equations 

(Ln(r’,s),  r(r,r',s))  = An(s)Ln(r,s)  (2.7a) 

and 

(?(?,?’, s),  Rn6',s))  = Xn(s)Rn(?,s)  , (2.7b) 


where  {^^(s)}  forms  the  set  of  eigenvalues  for  the  operator  defined  by 
f(s).  We  then  proceed  by  expanding  "everything  in  sight"  in  (2.2)  in 


terms  of  the  R^.  Then  by  means  of  the  biorthogonality  relation  between 
and  Rr,  we  arrive  at  the  explicit  representation 


- (l  (r'.s),  E(r',s)) 

J(r,s)  = l A (s)  —yz — r — — R (r,s)  . (2.8) 

n (Ln(r',s),  Rn(r',s))  n 


Thus,  we  may  write  the  inverse  kernel  T ^ in  (2.4)  as  an  expansion  of 
dyadic  terms 


T ^r.r'  ,s)  - £ X^U) 
n 


R (r , s)  L (r',s) 
n n 

(Ln(r',s),  Rn(r',s)) 


(2.9) 
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Now  the  poles  of  T ^(r,r',s)  are  observed  to  be  zeros  of  {X^Cs)}.  In 

[2]  Baum  discusses  the  specific  associations  between  elements  of  the 

set  of  poles  of  r ^ {s  } and  elements  of  the  set  of  eigenvalues  {X  (s)} 

m n 

and  designs  an  indexing  scheme  which  explicitly  represents  these 
associations.  We  are  not  concerned  with  the  specifics  of  the  association 
here,  however. 

Thence,  for  t such  that  the  asymptotic  behavior  of  J (r ,s)exp (st) 
allows  left  closure  and  assuming  simple  poles  for  notational  simplicity 


1 . _ _i  (l  (r',s),  E(r',s))  - _ 

J(r>t)  = j-7  f l * (s)  ~r= T Rn(r,s)  e ds 

2lfJ  Cfi  n n <Ln(?',s),  Rn(r',s))  n 


= l 

m 


dX 


As)  ~1  <L  (r'.s  ),  E(r',s  )) 
n m \ n m m r 


ds 


(l  (r'.s  ),  R (r',s  ))  n 
n n m n m / 


s t 

R (r,s  ) e m . (2.10) 
~ m 


The  summation  embraces  the  entire  pole  set  (s  }.  The  notation  used 

m 

assumes  that  X (s  ) = 0.  At  the  poles  the  eigenvectors  L (r,s  ) and 
n m n m 

R (r,s  ) become  the  coupling  vectors  C (r)  and  the  natural  mode  vector 
n ’ m r e>  m 

M (r) , respectively,  i.e.,  from  (2.7) 
m 


(C(r'),  *(?,?', s ))  =0 

\ m m / 


(2.11a) 


and 
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(r(r,?\s  ),  M (?'))  = 0 

> mm/ 


(2.11b) 


Using  this  modal  notation  and  the  relationship 

^n(Sm)  <Ln('*8m>*  to  •8m)  ’ Vr'  ’Sm^> 


dA 


ds 


(r'»8_),  R (r,s  )) 
' n m n m / 


we  rewrite  (2.10)  as 

\ __st 

J (r  jt)  = l B <C  (r),  E(r,s  ))  M (r)  e m 
m \ ra  m / m 

m 

The  normalization  constants 


6 * 


< 


c (r)»  r(r,r',s  ),  M (r')) 

m qs  mm/ 


-1 


(2.12) 


(2.13) 


are  introduced  in  arriving  at  (2.12).  The  inner  product  in  (2.12)  is 
generally  defined  as  a "coupling  coefficient"  depending  on  a particular 
spatial  form  of  E. 

The  formal  expansion  of  (2.12)  is  the  fundamental  SEM  form  - the 
singularity  expansion.  The  construction  by  which  we  arrived  at  (2.12) 
is  a formal  one.  As  subsequent  sections  of  this  chapter  show,  it  is 
indeeu  valid  for  a sufficiently  large  t.  For  t less  than  some  time  Tq , 
the  right-hand  closure  to  (2.6)  is  valid  and  J(r,t)  is  zero.  However, 
for  some  solutions,  there  is  a time  interval  during  which  neither  a zero 
value  nor  the  particular  expansion  (2.12)  is  a correct  representation  for 
J(r,t).  The  SEM  approach  for  these  solutions  seeks  to  modify  the  cou- 
pling coefficient  term  in  such  a way  that  a form  similar  to  (2.12)  holds 
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for  all  time  for  which  J(r,t)  is  nonzero.  In  succeeding  sections,  we 
present  some  alternative  coupling  coefficient  forms  to  this  end  and  basic 
guidelines  for  using  these  forms. 


2. 3 An  Analytically  Tractable  Example:  The  Transmission  Line 

2.3.1  Transmission  line  formulation 

The  lossless  transmission  line  proves  to  be  an  instructive  example 
in  studying  the  principle  involved  in  the  SEM  representation  . The 
instructiveness  of  the  problem  results  from  the  facts  that  the  SEM  rep- 
resentation is  fully  derivable  analytically  and  that  we  may  represent 
both  resonant  and  nonresonant  structures  with  finite  and  infinite  extent 
lines.  Because  the  residue  series  is  analytically  summable,  it  is  pos- 
sible to  explore  its  asymptotic  behavior  and  the  consequences  of  trun- 
cation of  the  representation  to  a finite  collection  of  poles.  It  happens 
that  the  problem  is  best  addressed  through  a differential  equation 
rather  than  an  integral  equation  of  the  form  (2.1).  However,  concepts 
applicable  to  the  Green's  function  here  transfer  directly  to  the  "inverse 
kernel"  Green's  function  of  the  integral  equation. 

The  geometry  of  the  transmission  line  is  defined  in  Figure  2.2. 

The  line  is  shorted  at  x = 0 and  x = L.  The  wave  equation  describing 
this  line  and  the  associated  boundary  conditions  are 


1_ 

C2  at2 


lU.t)  ■ - 

h uo-c)  -fe10-0  ■ °> 


I(x,t)  = 0, 


t < 0 


(2.14) 


The  author  is  grateful  to  Professor  D.  R.  Wilton  of  the  University  of 
Mississippi  for  his  suggesting  the  transmission  line  example  and  for 
providing  notes  on  parts  of  the  formulation  presented  here. 
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I 


Figure  2.2.  Defining  geometry  for  the  transmission  line  example. 
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In  the  transform  domain 


I(x,s) 


- ~ E(x,s) 

V 


(2.15) 


with 


h * <°-s>  -hl  a-s>  ■ 0 


We  may  obtain  a causal  Green's  function  solution  to  either  (2.14)  or 

(2.15) ,  and  the  two  are  a Laplace  transform  pair.  In  the  SEM  spirit  we 
write  an  s-domain  Green's  function,  i.e.,  the  Green's  function  for 

(2.15) .  This  function  is  defined  by 


~ (s/c)2  G(x,x',s)  = - 6 (x  - x')  , (2.16) 

ldxZ  J Z0C 


with 

G'(0,x',s)  = G'(L,x',s)  = 0 , Re  s > 0 

The  time  domain  counterpart  G(x,x',t  - t')  has  as  a boundary  condition 
G(x,x',t  - t')  » 0,  for  t < t'.  This  is  a causality  condition,  which 
appears  implicitly  in  (2.16)  by  way  of  the  region  of  validity  re  s > 0 
That  is,  we  must  have  a Green's  function  representation  which  is  valid 
on  the  Bromwich  contour,  the  real  part  of  which  is  taken  greater  than 


zero . 
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2.3.2  Solution  for  the  Green's  function 

We  choose  to  write  a constructive  solution  to  (2.16)  as 


G(x 


,x\S)  (r8lx-x’l  /c+l  e-«2nL/c 


^ l 


n=0 


-s[x+-x']/c  , -s[2L+(x-x? ) ]/c 

e + e + 


_>-s[2L- (x+x*  ) ]/c  + e~s  [2L-  (x-xf ) ] /c 


•])• 


Re  s > 0 . (2.17) 


This  form  is  taken  because  of  its  instructiveness  in  showing  explicitly 
the  traveling-wave  components  corresponding  to  a direct  wave  and  an 
infinite  summation  embracing  reflections  and  multiple  reflections  at  the 
shorted  ends  of  the  line. 

A closed-form  representation  results  directly  from  (2.17)  by  way  of 
the  sum  formula 


r -s2nL/c  ri  -2sL/c.  ^ _ 

l e = [ 1 - e ']  , Re  s > 0 

n=0 


The  closed-form  expression  is 


G(x,x’ ,s) 


1 / cosh  s [L  - |x-x’|]/c+  cosh  s [L  - (x  + x’ ) ] /c\  ia. 

2ZQ  ^ sinh  sL/c 


Alternatively,  this  form  can  be  constructed  directly  from  (2.16)  without 
using  the  summation  form.  The  expression  (2.18)  is  the  unique  analytic 
continuation  of  the  expression  (2.17). 

» — 1 

This  Green  s function  is  the  counterpart  to  the  inverse  kernel  T 
developed  for  the  integral  equation  formulation  in  Section  2.2.  For  this 
example,  however,  we  have  an  explicit  closed-form  representation  for  G. 

This  function  has  the  residue  series  representation 


G(x,x' ,s) 


Z L ^ en 
0 n*0 


s cos (mrx/L) cos (mix* /L) 
s^  + (nnc/L)^ 


(2.19) 
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where 


is  Neumann's  number.  It  may  be  shown  that  (2.19)  is  a complete  repre- 
sentation for  G(x,x',s)  as  expressed  in  (2.18).  That  is,  G(x,x',s) 
does  not  have  an  entire  function  component  in  its  Mittag-Lef f ler  repre- 
sentation. It  is  clear,  however,  from  (2.18)  that  asymptotically 
G(x,x',s)  is  exponential  in  character. 

We  also  observe  that  (2.19)  fits  the  dyadic  form  of  (2.9)  for  the  in- 
verse operator,  namely,  dyadic  numerator  terms  with  poles  in  the  denom- 
inator at  s • ±jmrc/L.  Thus,  for  a general  excitation  function  E(x',s), 
we  can  construct 


lo  e«< 

n*=0 


s cos (mrx/L) cos (mrx' /L) 
s^  + (nirc/L)^ 


2.3.3  The  causal  time  domain  Green's  function 
Alternatively,  we  can  invert  (2.19)  to  obtain  the  causal  time 
domain  Green's  function.  Then  a general  solution  may  be  formed  from 
the  convolution 


oo 

I (x,  t)  = / ^G(x,x',t  - t'),  E(x',t')^  dt' 

— 00 


(2.21) 


We  proceed  to  construct  G(x,x',t)  here  because  it  is  instructive  to 
observe  the  time  domain  effects  of  operations  which  we  perform  in  de- 
termining I(x,t)  from  the  transform  domain  form  (2.20). 
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The  time  domain  Green's  function  is  given  by 

G(x,x',t)  = / G(x,x',s)  est  ds  . (2.22) 


In  order  to  evaluate  this  integral  in  the  conventional  fashion,  we 
require  a knowledge  of  the  asymptotic  behavior  of  G.  This  is  discernible 
from  (2.18),  viz., 


G2(x,x' ,s) 


-s | x-x' | /c 
s|x-x' | /c 


in  r.h.p. 
in  l.h.p. 


(2.23) 


We  draw  attention  to  the  fact  that  this  information  is  obscured  in  the 
residue  series  form  (2.19).  It  is  also  noteworthy  that  G has  different 
asymptotic  forms  in  the  right-  and  left-half  planes.  Thus,  for 
t < |x  - x'|/c  we  may  close  Cg  to  the  right, and  since  G has  no  poles 
for  re  s > 0,  G(x,x',t)  = 0,  t < |x  - x'|/c.  For  t > -|x  - x’|/c  we 
may  close  Cg  to  the  left  and  obtain  the  residue  contributions  at  the 
poles  of  G(x,x',s).  We  write  then 


G (x, x' , t) 


< |x  - x'  | /c 

cos  (nirx'  /L)cos  (nirx/L)cos  (mrct/L)  , 


t > - | x-x' | /c  . 


(2.24) 

For  t e (-|x  _ x'|/c,  |x  - x'|/c)  the  function  must  be  zero,  but  the 
series  representation  is  valid,  too!  It  must  sum  to  zero.  Indeed,  the 
time  |x  - x'|/c  is  the  arrival  time  for  an  incident  wave  at  the  obser- 
vation point  x launched  from  a spatially  impulsive  generator  at  x' . 

The  phenomena  coming  to  play  are  seen  if  we  use  the  following  distri- 
butional identity  to  sum  the  above  series. 
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£ 6 (x  - an)  - £ eJ2nnx/a  «,  £ e cos(2nirx/a) 

_oo  _oo  n-0  n 

This  yields 

oo 

G(x,x',t)  = ~^z  £ ( 6 [ t - (x  - x')/c  - 2nL/c] 

0 n=-»> 

+ 6[t  + (x  - x')/c  - 2nL/c] 

+ 6[t  - (x  + x')/c  - 2nL/c] 

+ 6 [ t + (x  + x')/c  - 2nL/c]},  t > -|x  - x'j/c  , 

(2.25a) 

or 


G(x,x' ,t) 


I x - x 1 | / c ] + 6 [ t - (x  + x')/c] 


+ £ {6[t  - |x  - x'|/c  - 2nL/c] 

n-1 

+ 6 [ t + |x  - x'|/c  - 2nL/c] 

+ 6[t  - (x  + x')/c  - 2nL/c] 

+ 6 [ t + (x  + x')/c  - 2nL/c)}|  . (2.25b) 

This  distribution  is  represented  symbolically  as  a function  of 
time  in  Figure  2.3.  The  distribution  represents  a direct  current  impulse 
arriving  at  the  observation  point  at  a time  |x  - x'|/c  accompanied  by 
three  impulses  arriving  after  undergoing  reflections  at  one  or  both 
shorted  ends  of  the  line.  This  pattern  is  replicated  with  a period 
2L/c  since  the  line  is  lossless  and  nondispersive . 


We  observe  that  the  series  representation  yields  replication  into 
negative  time  as  well,  as  shown  by  the  dashed  line  impulses  in  Figure  2.3. 
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It  is  clear  that  if  we  attempted  to  apply  the  series  for  t < -|x  - x'|/c, 
we  would  obtain  erroneous  results.  The  phenomenon  of  having  both  the 
series  representation  and  the  zero  representation  for 

-|x  - x* |/c  < t < |x  - x'|/c  is  evident,  too.  The  series  converges  to  a 
distribution  which  has  zero  weight  on  this  interval. 

A "physical"  interpretation  may  be  placed  upon  the  acausal  pulses 
which  appear  in  (2.25a).  Figure  2.4  is  a diagram  showing  the  spatial 
location  of  the  current  impulses  represented  by  (2.25a)  for  a source 
generator  6 (x  - x')6(t).  The  locations  are  shown  for  four  times  - spaced 
at  an  increment  At  - with  the  first  sample  taken  at  -At.  Using  the 
t = 0 (second)  figure  as  a reference,  we  see  that  the  negative  time 
impulses  in  the  first  figure  are  simply  impulses  propagating  toward 
the  source  so  as  to  coalesce  at  the  source  point  at  t = 0.  All  of  the 
other  negative  time  terms  are  earlier  reflections  of  these  pulses.  The 
subsequent  propagation  away  from  x = x'  and  reflection  at  the  short  at 
x = L shown  in  the  third  and  fourth  figures  are  evident.  The  current 
impulses  occurring  in  negative  time  are  a result  of  applying  the  singu- 
larity expansion  for  G(x,x',t)  too  early. 

2.4  An  Examination  of  the  Time  Domain  Green's  Function  for  the  Thin-Wire 

Scatterer 

The  preceding  section  describes  the  analytical  aspects  of  an  SEM 
solution  for  a shorted  section  of  lossless  transmission  line.  The  con- 
sequences of  the  analysis  related  to  the  inversion  are  evident  when 
the  solution  is  viewed  in  terms  of  the  time  domain  Green's  function 
for  the  problem.  We  now  turn  to  the  scattering  problem  perhaps  most 
closely  associated  with  the  transmission  line  - that  of  the  scattering 
from  a perfectly  conducting  thin  wire. 


x*  0 

Figure  2.4. 


Diagram  of  an  acausal  Green's  function  for  the  shorted 
transmission  line.  The  spatial  location  of  the  impulses 
of  current  on  the  line  are  shown  "frozen"  in  time  for 
four  times  each  advanced  At  over  its  predecessor.  Time 
t - 0 is  given  in  the  second  graph. 
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The  SEM  solution  to  this  problem  must  be  constructed  numerically, 
however.  As  a consequence,  we  have  available  only  a limited  number  of 
terms  in  the  singularity  expansion.  There  is  no  analytical  basis  from 
which  we  can  observe  the  asymptotic  behavior  of  the  transform  domain 
Green's  function  (inverse  kernel)  for  this  problem.  We  can,  however, 
compute  it  within  the  limits  of  approximation  dictated  by  the  amount  of 
available  pole  information.  Through  the  previous  example  we  attempt 
to  conceptually  relate  features  of  the  time  domain  and  asymptotic  transform 
domain  behaviors.  Practically,  we  determine  how  to  construct  the  singu- 
larity expansion  from  directly  discernible  time  domain  ramifications. 

We  rely  on  the  results  of  Tesche  [8], [9]  for  numerically  determined 
SEM  quantities  for  the  thin  wire.  A tabulation  of  poles,  coupling  co- 
efficients, and  normalization  constants  along  with  graphs  for  the  natural 
mode  currents  are  given  in  [9]  for  ten  so-called  first-layer  poles  and 
three  second-layer  poles.  The  cylinder  for  which  these  data  are  pro- 
vided has  a diameter-to-length  ratio  of  0.05. 

The  geometry  under  consideration  is  defined  in  Figure  2.5.  The 
wire  is  located  along  the  x-axis  extending  from  -h  to  h.  We  construct 
the  inverse  kernel  for  the  structure  from  (2.12)  by  taking 
E(r,t)  * u^  6 (x  - x' )6 (t  - t ' ).  By  using  the  usual  thin-wire  analysis 
concepts,  we  observe  that  the  coupling  vectors  and  the  currents  are 
scalar  functions  of  a scalar  variable  x.  Further,  each  coupling  vector 
is  equal  to  its  associated  current  mode.  Accordingly,  we  write 

-i  / -s  t'  s t 

T i(x,x',t  - t')  = u(t  - t')  T B (J  (x),  S(x  - x')  e m ) J (x)  e m 

m 

s (t-t') 

- u(t  - t')  l em  J (x’)  J (x)  m 
u m m m 

m 


e 


(2.26) 
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We  include  the  unit  step  function  u(t  - t')  at  this  stage  as  a pre- 
caution to  insure  that  our  inverse  kernel  is  causal.  This  matter  is 
explored  further  below.  Then  for  an  arbitrary  incident  field,  we  con- 
struct the  current  response  as 

oo 

J(x,t)  = / (r  ^(x.x'.t  - t'),  u^  • E(x',t')^  dt'  . (2.27) 

— CO 

The  motivation  for  exploring  this  example  is  to  show  that,  indeed, 
an  acausal  response  is  contained  in  the  expansion  in  (2.28)  and  the  unit 
step  function  is  necessary  to  enforce  causality  in  the  inverse  kernel. 

This  is  readily  tested  by  expanding  F ^ for  specimen  source  coordinates 
x',  t'  and  studying  the  behavior  of  this  expansion  for  positive  and 
negative  times.  A plot  of  a test  expansion  of  T ^ for  the  wire  is  given 
in  Figure  2,6.  The  expansion  is  taken  with  the  source  point  x'  =0, 
t'  = 0 for  two  different  observation  points  x = 0 and  x = 0.6  h. 

We  observe  from  the  figure  an  impulsive  current  at  x = 0 and  t = 0 
due  to  the  generator’s  "firing."  This  impulse  divides  and  propagates 
outwardly  to  produce  impulses  at  ±0.6  h beginning  at  a time  0.6  h/c. 

The  reflections  of  these  pulses  are  seen  at  ±0.6  h,  when  t = 1.4  h/c, 
and  then  merging  at  2 = 0 when  t = 2 h/c.  This  qualitatively  predictable 
reflection  pattern  continues.  We  observe  an  oscillatory  current  function 
between  the  expected  pulses.  This  behavior  is  due  to  the  truncation 
of  the  residue  expansion  to  a finite  number  of  poles.  In  fact,  as 
mentioned  in  Section  2.2.1,  T ^ cannot  be  a function,  it  must  be  a 
Dirac-like  distribution.  The  smoothness  and  the  oscillations  are 
attributable  to  expansion  of  the  distribution  in  terras  of  a finite 


number  of  continuous  functions. 


L»WV 


Figure  2.6.  Time  domain  inverse  k 
as  a function  of  time 
observation  points, 
data  at  ten  first  lay 


t { h/c) 


ernel  for  thin-wire  scatterer  plotted 
for  a source  point  at  x'  * 0 and  two 
The  function  was  expanded  from  SEM 
er  poles  as  given  by  Tesche  [9]. 
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Ferhaps  the  most  important  feature  to  observe,  however,  is  the  pres- 
ence of  an  impulse  at  x = 0,  at  t = -2  h/c  and  at  x = ±0.6  h,  t = -0.6  h/c. 
These  impulses  are  acausal!  We  see  that,  indeed,  we  must  impose  causality 
on  T * with  a gating  function  in  (2.26).  The  residue  series  form  for 
the  inverse  kernel  contains  impulses  in  negative  time  which  are  acausal. 

The  effects  of  these  pulses  in  a finite  series  become  more  pronounced  with 
increasingly  negative  time  as  one  would  expect  from  terms  of  the  form 
exp(s^t)  with  re  s^  < 0.  Other  data  not  shown  herein  indicate  that  the 
inclusion  of  second  layer  poles,  whose  real  parts  are  more  negative,  causes 
a very  rapid  blowup  in  negative  time. 

The  physical  interpretation  of  the  negative  time  impulses  is  es- 
sentially the  same  as  that  given  for  the  transmission  line  in  the  preceding 
section  if  we  consider  the  additional  effects  of  dispersion  and  radiative 
loss  in  the  antenna  problem. 

The  negative  time  impulses  in  the  SEM  representation  for  the  time 
domain  Green's  function  affects  the  form  of  the  coupling  coefficient  in 
the  early  time  of  the  transient  response  of  a scatterer.  We  discuss 
this  bearing  in  Section  2.5.3. 

2 . 5 Relation  Between  Coupling  Coefficient  Forms  and  Asymptotic  Behavior 
2.5.1  Effect  of  closure  time  on  the  coupling  coefficient 
In  Section  2.2.3  it  is  indicated  that  the  so-called  coupling  coef- 
ficient for  the  singularity  expansion  results  when  t is  sufficiently 
large  that  we  can  close  the  contour  to  the  left  in  the  Laplace  inversion. 

We  now  turn  to  the  details  of  the  interrelation  between  the  asymptotic 
form  for  the  transform  domain  current  solution  and  specific  forms  of 
the  coupling  coefficient. 
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We  begin  with  (2.6),  the  expression  for  the  current  solution  to 
an  electromagnetic  integral  equation. 

( ?-1  (r,r ' ,s) , E(r,s))  eSt  ds  . (2.28) 

CB 

For  the  present  discussion,  we  consider  the  incident  field  E(r,t)  to 
be  a plane  wave  with  a time  history  f(t).  We  restrict  f(t)  to  have  a 
transform  f(s)  which  is  algebraic.  Further,  we  consider  object  poles 
alone.  The  concepts  which  follow  are  readily  extended  to  include  wave- 
form poles  and  more  general  f(s).  These  restrictions  are  cast  in  the 
interest  of  notational  simplicity.  A plane  wave  field  may  be  written 
in  the  time  domain  as 

E (r , t)  = Eq  f(t)  * 6 ( t - p • r/c)  , (2.29) 

where  Eq  is  the  polarization  vector  of  the  wave,  (5  is  a unit  vector  in 
the  direction  of  wave  propagation  and  denotes  convolution.  This 
representation  has  the  transform  domain  image 

E(r,s)  = Eq  f(s)  e r^c  . (2.30) 

We  use  this  expression  in  (2.6),  change  the  order  of  integration  between 
the  spatial  and  transform  variables,  and  denote  the  spatial  integration 
explicitly  to  obtain 

J(r,t)  “ / / ” 1(r,r’,s)  • E„  e r f(s)  eSt  dsdr'  . (2.31) 

3 Body  CB 
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The  interchange  of  integration  operators  allows  us  to  observe  the 
asymptotic  behavior  of  the  integrand.  Were  the  spatial  integral  per- 
formed first,  the  asymptotic  behavior  would  be  obscured  in  the  process. 

The  s-dependent  terms  in  the  integrand  in  (2.31)  are  the  inverse 
kernel,  the  exponential  of  the  propagating  incident  wave,  the  trans- 
form of  the  incident  wave  time  history  f(s),  and  the  Laplace  kernel. 

as—  X 

The  asymptotic  behavior  of  T is  generally  unavailable  to  us.  The 
propagating  wave  term  and  the  Laplace  kernel  each  have  explicit  ex- 
ponential behavior.  Thus,  the  algebraic  behavior  of  f(s)  and  the 

5-  X 

potential  algebraic  behavior  of  T are  dominated  by  the  exponential 
factors.  This  restricts  our  interest  in  the  asymptotic  behavior  of 

F * to  exponential  terms.  Accordingly,  we  define  asymptotic  forms  for 

+ . " 

C and  C as 

QQ  <30 


T_1(r,r' ,s)  - 


fm  -sTR(?,r')  + 

Dr  dR(s)  e s e 


LDL  dL(s)  6 


-sTL(r,r') 


(2.32) 


e C 


where  the  d(s)  functions  are  algebraic  and  the  D are  dyads  wtiich  are 
independent  of  s.  The  coefficients  of  s in  these  terms  are  dimensionally 
equal  to  time.  They  depend  on  both  r and  r',  potentially,  and  might  be 
positive  or  negative. 

We  restrict  attention  to  cases  where  T,  (r,r')  < T (r,r')  so  that 

L K 

pole  information  suffices  in  the  representation  for  the  Laplace  in- 
version of  r 
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For  t > T^(r,r')  the  Integrand  in  (2.31)  decays  exponentially  as 

Is  I 00  on  Coo*  the  left-side  closure  of  CD.  For  such  t we  close  C 

D B 

with  and  a residue  series  results  from  the  inner  integral.  For 
c < ^L^r,r  ) — ^ r * r ' ) the  integrand  decays  on  C+  so  that  closing  to 

the  right  yields  zero.  In  the  dyadic  notation  of  (2.12),  we  write 
explicitly 


J(r,t) 


/ u[t  - p • r'/c  - T (r,r' ) ] \ 
Body  m=0 


B E 
m C 


_ _ _ S (t-^T'/c) 

C (r')f(s  ) e m 
m m 


* Mm(r)  dr'  , (2.33) 


where 


u(t)  = 


A bandlimited  behavior  in  f(s)  allows  us  to  truncate  the  series  at, 
say,  M terms  in  (2.33).  Then  with  some  rearrangement  we  have 


3(?.t)  - J f(..)  6.  / I • C S') 

m Body 


-s  p-r'/c 


u[t  - B * r'  - TL(r,r')]  dr* 


M (r)  e 
m 


s t 
m 


(2.34) 


The  integral  factor  is  a "generalized"  coupling  coefficient 


n”rr,t)  ' Body  ’E°  ' E"rr,) 


-s  e-r'/c 


u[t  - p • r'  - TL(r,r')]  dr' 


(2.35) 
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Then  (2.34)  is  written  as  the  singularity  expansion  form 

M s t 

J(r,t)  - 7 f(s  ) 6 n (^.t)  M (?)  e m 
**  m m m m 

m 


(2.36) 


There  are  two  commonly  u ed  coupling  coefficient  forms  which  may  be 
interpreted  as  parti  ulari/at i ns  of  the  generalized  form  in  (2.35). 

Baum  terms  these  as  cl ■ 1 and  class  2 coupling  coefficients,  respectively 

[10]. 

A class  1 form  results  by  assuming 

TL (?,?’)  - -p  • r'/c  + t'  (2.37) 

where  t'  is  some  chosen  "turn-on  time"  typically  taken  to  be  the  arrival 
time  of  the  incident  wave  at  the  body.  This  results  in  the  current 
expansion 

M _ n x st 

J(r,t)  = u(t  - t')  l i(s  ) 8 M (r)  e m , (2.38) 

u m m m m 
m 


where  the  class  1 coupling  coefficient  is 


Body 


m 


-s  § 
m 


r'/c 


dr' 


(2.39) 


It  is  notable  that  the  assumed  T form  is  independent  of  r and  yields  a 

la 

coupling  coefficient  independent  of  both  r and  t.  For  some  applications 
of  SEM  such  as  broadband  equivalent  circuit  modeling  of  distributed 
structures  [11],  this  feature  is  useful. 

A class  2 coupling  coefficient  results  when  T is  taken  to  be  zero. 

la 

Then  a series  in  the  form  (2.36)  results  but  with  the  coupling  coefficient 
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-s  p*r'/c 

/ E • C (r')  e m u(t  - p • r'/c)  dr' 

„ , U m 
Body 


* / 

R(t) 


C (r') 
m 


-s  f>*r'/c 
m 
e 


dr' 


(2.40) 


The  integral  over  R(t),  the  "excited  region,"  is  explained  with  the  aid  of 
Figure  2.7.  The  effect  of  the  unit  step  u(t  - p • r'/c)  in  the  coupling  co- 
efficient integral  is  to  gate  the  integrand  so  that  it  contributes  to  n (t) 

n 

only  over  the  part  of  the  body  at  which  the  incident  field  has  arrived.  A 
notable  feature  of  the  class  2 form  is  that  it  is  time  dependent.  After 
the  leading  edge  of  the  incident  wave  has  cleared  the  object,  it  takes  on 
a constant  value  identical  with  the  class  1 value. 

To  the  best  of  our  knowledge,  previous  workers  in  SEM  have  made  their 
choice  of  coupling  coefficient  forms  without  benefit  of  the  asymptotic  form 
of  the  inverse  kernel.  In  fact,  the  class  2 choice  results  from  an  assumption 
that  the  inverse  kernel  has  no  exponential  asymptotic  behavior,  i.e.,  T = 0 

Li 

as  in  the  above.  (An  equivalent  point  of  view  is  termwise  integration  of  the 
series.)  The  class  1 coefficient  is  a choice  of  convenience  in  the  form  of 
time  independent  coupling  coefficients.  In  the  subsections  which  follow,  we 
explore  some  ways  of  reconciling  coupling  coefficient  forms  with  their  inverse 
kernel  expansions. 

2.5.2  Heuristic  development  of  the  asymptotic  form  of  T - a 
conservative  coupling  coefficient 

The  transmission  line  example  developed  in  Section  2.3  gives  us 
insight  because  the  asymptotic  form  of  the  transform  domain  Green's 
function  is  clearly  discernible.  In  that  example,  there  is  some  freedom 
of  choice  in  the  left  closure  time  because  of  the  time  interval  during 
which  both  the  SEM  representation  and  zero  are  valid  representations  of 


the  solution. 


Figure  2. 
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Leading  edge 
of  incident  wave 


. Representation  of  a scatterer  showing  the  "excited  region" 
over  which  the  illuminating  wave  has  passed. 
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The  other  example  in  our  knowledge  where  the  asymptotic  form  of 
the  SEM  Green's  function  has  been  explored  is  in  the  loop  [5].  For 
Lue  majority  of  SEM  problems,  which  rely  on  numerical  determination  of 
the  component  terms  for  the  expansion,  it  seems  unlikely  that  the 
necessary  asymptotic  information  will  be  directly  available.  In  the 
following  paragraphs,  we  develop  heuristically  a useable  form  for  the 
asymptotic  behavior  and  a new  coupling  coefficient  form  which  results 
from  it.  It  is  observed  to  be  unnecessarily  conservative  for  the 
examples  of  the  transmission  line  and  the  wire  scatterer  considered 
in  previous  sections.  But,  because  of  its  conservatism,  it  provides  a 
"safe"  approach  to  new  problems.  The  next  subsection  provides  a nu- 
merical test  procedure  for  relaxing  the  conservatism  to  an  appropriate 
degree . 

We  consider  the  asymptotic  behavior  of  the  inverse  kernel  with 
the  aid  of  Figure  2.8.  The  object  is  represented  in  a coordinate  system 
with  origin  at  0.  A source  point  is  located  at  r = r'  and  an  observa- 
tion point  at  r « r^.  The  vector  s is  a local  surface  vector  at  r' 
with  unit  components  in  two  orthogonal  directions  in  the  local  tangent 
plane.  At  the  observation  point  r^,  the  current  on  the  body  is  repre- 
sented by  J(rg,t)  in  local  vector  components.  Then  if  an  excitation 

E(r,t)  - S(r')6(r  - r')6(t)  , 

which  is  impulsive  in  time  and  space,  is  applied  at  the  observation 
point,  we  expect  a temporal  response  which  departs  from  zero  only  for 
t > | r — r * | / 


c.  That  is 
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J(rQ,t) 


I (rQ,t) 


c < |r0  - *' |/c 
t > 1*0  “ r'l/c 


with  I(rQ,t)  = 0 only  at  isolated  points  in  t. 

By  definition,  r ^(r,r’,t)  is  the  kernel  of  a convolution  relating 
E(r,t)  and  J(r,t)  so  that 

oo 

j(r0,t)  = / (r"1(?0,?,t  - t'),  S(r ' )6 (r  - r')6( t'))  dt' 

—00 

= ^_1(?0,^’,t)  • s (r » ))  . 

It  follows  that  in  the  transform  domain  J(Tq,s)  and  hence  T-1 (r^ , r ' ,s) 
must  have  an  exponential  factor  to  represent  the  time  delay  | r^  - r'|/c. 
We  write  explicitly 


5-1 ~ i N -s  r-r ' 
r (r,r  ,s)  = e 1 


/C  KU.'r'.s) 


and 


H(?,?\s)  - D e-Sl?-r'l/c 


(2.41) 


where  K is  a dyadic  function  and  5 is  a dyad  independent  of  s. 

This  assumed  asymptotic  behavior  implies  a closure  time  term  of 
the  form 


TL(r,r)  = - | r - r' | /c 


in  (2. 32)- (2. 35) . Specifically,  a coupling  coefficient  of  the  following 
form  results. 
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in 


-s  p*r'/c 

/ Eq  • Cm(?')  e m u[t  - 0 • ?' 

Body 


|r  - r* | /cl  dr'  . 


(2.42) 

We  term  this  a "class  X"  coupling  coefficient.  (An  X-tra  cautious  one.) 

Some  observations  are  in  order.  First,  the  class  X coupling  co- 
efficient is  dependent  on  both  the  observation  time  and  the  observation 
point.  This  feature  renders  it  more  complex  than  either  the  class  1 or 
class  2 forms  discussed  in  Section  2.5.1.  Second,  the  effect  of  the  unit 
step  function  in  (2.42)  is  to  enforce  causality  directly.  This  feature  is 
seen  clearly  in  Figure  2.9  for  the  example  of  a thin-wire  scatterer.  The 
figure  represents  the  support  of  current  induced  by  the  incident  wave  at 
two  coupling  points:  x = L/4  and  x = 3L/4.  The  coupling  is  from  an 

incident  plane  wave  whose  direction  of  propagation  forms  a"  angle  of 
60  degrees  with  the  axis  of  the  scatterer.  Consequently,  its  velocity 
component  along  the  scatterer  is  2c.  The  time  sequence  depicts  the 
wave  arriving  at  the  scatterer  at  t = 0,  coupling  at  x = L/4  t = L/8, 
progressing  on  at  t = L/4c,  coupling  at  x = 3L/4  t = 3L/8,  etc.  The 
support  for  the  current  induced  when  the  wave  is  at  x = L/4  is  depicted 
as  impulsive  when  the  wave  is  at  that  point  then  propagating  outwardly 
at  the  velocity  of  light.  A similar  phenomenon  is  pictured  for  coupling 
at  x * 3L/4. 

A third  observation  is  that  this  explicit  enforcement  of  causality 
may  prove  unduly  stringent.  For  example,  the  resonant  transmission 
line  problem  in  Section  2.3  led  to  a Green's  function  with  the  asymp- 
totic behavior 


C(x,x',s)  - 


in  r.h.p. 
, in  l.h.p. 


incident 
Plane  wave 

Ee8(t-rp/c) 


1 -*  1 - I I 

t=0  X=U/4  X = U/2  X»  3L/  4 X=L 


^ ♦ * 

, > , 

1 

1 t a 0 

r yrV  1 

1 1 ?////{///, 

1 

1 t = L/8C 

1 4 • 1 / a r* 

Y/////////////7Z\ 

j T ' L/4  C 

1 t-  3 L/ec 

t * L /2C 


Figure  2.9.  Representation  of  the  expanding  of  support  of  current 

propagating  away  from  source  coupling  at  L/4  and  3L/4  as 
a plane  wave  passes  across  a linear  scatterer. 
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i.e.,  different  asymptotic  forms  on  and  so  that  there  is  a period 
of  time  during  which  either  closure  of  C,.  is  valid.  For  this  example 
we  can  choose  a left  closure  time  T'  = 0 and  obtain  a class  2 coupling 

L 

coefficient  as  in  (2.40). 

In  the  following  subsection,  we  relate  the  function  of  the  asymp- 
totic closure  time  to  the  time  domain  inverse  kernel.  Thereby,  we  de- 
termine tests  which  may  be  applied  to  the  numerically  expanded  time 
domain  Green's  function  to  determine  what  freedom  is  available  in  a 
particular  problem  in  choosing  a closure  time  and  concomitantly  the 
coupling  coefficient  form. 

2.5.3  The  purpose  of  the  closure  time  in  the  time  domain  Green's 
f unct ion 

In  this  section  we  present  some  numerical  tests  which  may  be  applied 
for  any  particular  form  to  determine  the  acceptability  of  a given  coupling 
coefficient  form.  We  develop  these  tests  based  primarily  on  insight  into 
the  relationship  between  closure  time  and  the  time  domain  form  of  the 
inverse  kernel  as  expanded  in  a singularity  expansion  which  we  gain 
from  the  transmission  line  and  straight-wire  scatterer  examples. 

To  begin,  we  recognize  the  equivalence  between  the  algebraic  multi- 
plication of  the  transform  domain  inverse  kernel  with  the  incident  field 
and  the  time  convolution  of  the  counterpart  time  domain  functions. 

That  is,  the  transform  domain  form  for  J(r,t)  of  (2.6) 

J(r,t)  * -JTT  / (f  '(r.r'.s),  E(r',s))  est  ds 
CB 


is  equivalent  to  the  time  domain  convolution 
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3(r,t)  = / (r_1(r,r,,t  - t'),  E(^'.t'))  dt'  . (2.43) 

— OO 

Any  assumed  asymptotic  behavior  applied  in  the  inversion  of  (2.6)  is 
carried  in  (2.43)  in  the  individual  inversions  of  F ^ and  E. 

To  illustrate  we  again  call  upon  the  straight-wire  scatterer  ex- 
ample. Suppose  that  we  expanded  its  response  due  to  an  obliquely  in- 
cident plane  wave  using  a class  1 coupling  coefficient.  In  Section  2.5.1 
Equation  (2.37),  it  is  pointed  out  that  the  class  1 coupling  coefficient 
assumes  an  asymptotic  form 

?-!/-  f.  -s(t'-p*r'/c) 

I (r,  r ,s)  ~ D e r , 

where  D is  a dyadic  constant.  Let  us  consider  a configuration  like 
that  in  Figure  2.9  and  take  the  turn-on  time  t'  = 0,  the  time  at  which 
the  incident  wave  first  impinges  on  the  wire.  The  class  1 assumption  is 


The  extreme  value  of  r'  in  the  spatial  inner  product  is  |r'|  = L.  In 
this  extreme  by  using  the  class  1 coupling  coefficient,  we  effectively 

s—  — — 

gate  the  singularity  expansion  of  T (r,r',t)  with  the  unit  step 
u(t  + p • r'/c)  with  |r'|  = L.  This  requires  that  the  singularity  ex- 
pansion manifest  causal  behavior  for  all  t > -f>  • r'/c.  In  the  60  degree 
incidence  example  this  translates  to  t > -L/2c. 

A reference  to  the  plotted  singularity  expansion  of  T ^(r,r',t) 
for  the  wire  given  in  Figure  2.6  indicates  that  the  required  behavior 
is  not  present.  There  is  an  acausal  impulse  at  -|x  - x'|/c  so  that 
the  required  clear  time  is  not  present  for  all  observation  paints  on 
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the  wire.  On  the  other  hand,  the  class  2 coupling  coefficient,  which 
5—1 

assumes  that  T has  no  asymptotic  exponential  behavior,  i.e.,  it  is 
algebraic,  dictates  a required  causal  behavior  in  the  expansion  for 
t > 0.  We  observe  from  the  figure  that  this  is  the  case  within  the 
frequency  and  numerical  error  limitations  of  the  expansion. 

Tesche  [12]  has  attempted  to  use  the  class  1 coupling  coefficient 
without  success  on  the  thin  wire.  The  results  which  he  presents  in 
[8], [9]  indicate  that  the  class  2 coefficient  form  is  satisfactory. 

Wilton  and  Umashankar  have  experimented  with  both  forms  on  the  L-wire 
scatterer  as  in  [5].  Their  experience  has  been  similar  to  Tesche's 
with  regard  to  the  success  and  failure  of  the  forms  [13].  This  experience 
is  consistent  with  the  theoretical  evidence  herein. 

This  illustration  suggests  a test  which  may  be  applied  in  the  time 
domain  to  determine  the  appropriateness  of  a given  coupling  coefficient 
form  to  a particular  problem.  Corresponding  to  any  coupling  coefficient 
form  is  an  assumed  asymptotic  behavior  for  the  transform  domain  inverse 
kernel.  This  asymptotic  form  in  turn  implies  a turn-on  time  for  the  time 
domain  inverse  kernel.  For  specimen  source  and  observation  points,  the 
time  domain  inverse  kernel  can  be  numerically  expanded  in  the  manner  of 
Section  2.4.  By  observation,  one  can  determine  whether  or  not  the  turn-on 
time  associated  with  the  coupling  coefficient  in  question  and  the  observed 
time  domain  inverse  kernel  are  consistent  with  causality.  Stated  as  a 
recipe  : 

A.  Analytical  steps 

1.  For  the  particular  coupling  coefficient  to  be  tested,  d*  ‘rn 
the  support  of  the  integrand  in  the  coupling  coefficient  inner 
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product,  e.g.,  u(t  - t')  as  in  (2.38)-(2.39)  for  class  1, 

u(t  - £ • r'/c)  in  (2.40)  for  class  2, 

u(t  - p • r'  - |r  - r’|/c)  in  (2.42)  for  class  X. 

2.  Translate  the  unit  step  argument  into  the  left-half  plane 
asymptotic  form  which  it  assumes,  i.e.,  exp(sA),  where  A is  the 
argument  of  the  unit  step.  For  class  2 plane  wave  excitation 
this  is 

s(t-p*r'/c)  st  -sp*r'/c  , 0. 

e 1 = e e (*e  ) 

St 

3.  Factor  out  the  e Laplace  kernel  and  the  incident  wave 
asymptotic  behavior.  For  class  2 this  takes  care  of  both  factors 
above.  The  remaining  term  - e^  above  - is  the  assumed  asymptotic 
behavior  of  the  inverse  kernel  implicit  in  the  coupling  coefficient 
form.  The  exponent  of  this  exponential  is  -s  times  the  implicit 
turn-on  time  of  the  time  domain  inverse  kernel. 

B.  Numerical  steps  (assuming  that  poles,  modes,  etc.,  are  known) 

1.  Pick  some  specimen  observation/source  point  pairs  representing 
extreme  cases  (e.g.,  maximum  separation  on  the  object)  and  the 
transition  in  between. 

2.  For  sources  impulsive  in  time  and  space,  expand  the  time  domain 
inverse  kernel  in  the  manner  of  Section  2.4. 

3.  Observe  the  "clear  time"  in  the  time  domain  kernel  prior  to 
the  arrival  of  the  first  causal  disturbance.  In  particular, 
note  the  dependence  of  this  time  on  r and  r' . 

C.  Compare  the  required  clear  time  with  the  turn-on  time  determined  in 

Step  A3.  Does  the  inverse  kernel  provide  clear  time  between  the  turn- 
on time  implied  coupling  coefficient  and  the  first  causal  disturbance? 
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If  the  answer  to  this  question  is  "No,"  then  a more  conservative  coupling 
coefficient  must  be  used.  The  class  X coefficient  developed  in  the 
preceding  subsection  is  highly  conservative  since  it  imposes  causality 
directly.  The  class  2 coefficient  is  somewhat  less  conservative  be- 
cause it  requires  a clear  time  from  t - 0 to  the  first  causal  disturbance. 
Experience  in  using  it  on  many  problems  corroborates  its  usefulness, 
however.  The  class  1 coefficient  should  be  used  only  with  caution. 

It  assumes  a clear  time  of  -L/c,  independent  of  r,r',  where  L is  the 
maximum  dimension  of  the  object. 

A second  potential  test  procedure  on  a coupling  coefficient  form  is 
simply  to  compare  its  results  with  the  conservative  class  X form.  One 
must  be  cautious  to  include  representative  incident  angles  and  spectra 
of  the  excitation  in  conducting  this  test,  however. 

The  Step  B3  above  is  the  point  at  which  one  must  exercise  proverbial 
"engineering  judgment."  Obviously,  in  dealing  with  graphical  data  such 
as  that  in  Figure  2.6,  the  decisions  are  not  clear-cut.  We  use  the 
qualitative  plot  of  a time  domain  inverse  kernel  in  Figure  2.10  to  point 
out  the  features  which  one  might  expect  in  such  a plot.  Hopefully, 
the  most  prominent  feature  will  be  the  expected  causal  disturbance  ar- 
riving at  t = |r  - r'j/c.  The  clear  time  region,  if  present,  must  be 
an  approximation  to  a zero-weight  distribution.  In  Figure  2.10  we  see 
depicted  oscillations  which  are  significant  in  magnitude  compared  with 
the  proper  disturbance.  We  note,  however,  that  the  waveform  with  which 
this  function  is  convolved  varies  slowly  compared  with  these  oscil- 
lations if  it  is  bandlimited  appropriately  for  the  number  of  poles  used. 
Thus,  in  an  approximate  sense,  these  ripples  cancel  among  themselves 
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rW.r'.t) 


Figure  2.10.  Qualitative  representation  of  a time  domain  inverse 

kernel  plot  showing  the  principal  features  which  have 
bearing  on  the  applicability  of  the  kernel. 


46 


in  the  convolution.  The  acausal  region  of  the  distribution  may  be  marked 
by  one  of  two  features:  either  a prominent  disturbance  prior  to 

t = |r  - r'|/c  or  a "blowup"  of  the  expansion.  This  blowup  is  expected 
at  some  point  in  negative  time,  because  in  negative  time,  the  real  parts 
of  the  exponents  in  the  series  are  positive  (the  result  of  the  negative 
real  part  of  the  pole  times  a negative  time). 

Another  breakup  of  the  representation  to  which  one  should  be  sen- 
sitive in  studying  F ^ is  that  which  is  expected  from  attempting  to 
apply  the  SEM  method  at  too  high  frequencies.  That  such  a breakup  is 
expected  is  based  on  the  following  reasoning.  In  the  dyadic  representa- 
tion for  the  time  domain  inverse  kernel  on  the  wire  scatterer,  for  ex- 
ample, 

, s t 

T (x,x',t)  = [ 8 e B J (x ) J (x') 
ro  m m 

To  a first  approximation,  the  function  of  the  spatially  dependent  factor 

J (x)  J (x')  is  to  delay  a pulse  in  time  by  a time  lx  - x'|/c,  viz., 
mm 

if  at  the  observation  point  x'  = x 
M s t 

T i(x,x,t)  = £ A e m = f(t)  , 

m*0  m 

at  a point  separated  from  the  observation  point 

M -s  j x-x' J / c st 

T (x,x',t)  > f(t  - |x  - x'|/c)  *=  £ A e m e m 

ra=0  m 

It  is  seen  that  if  M is  increased,  the  factors  exp[-sm|x  - x'|/c]  grow 
since  Re  s becomes  more  negative  as  m increases.  Further,  the  larger 
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a pole  in  magnitude,  the  more  subject  it  is  to  numerical  error.  This 
error  is  magnified  due  to  the  exponential  dependence  on  the  s^.  This 
error  will  manifest  itself  as  a blowup  in  the  oscillations  in  the 
clear-time  region  of  the  inverse  kernel.  The  net  effect  is  a numerical 
limitation  on  the  time  resolution  relative  to  body  size  in  a given  SEM 
solution. 

2 . 6 Conclusions 

In  this  chapter  we  present  the  fundamental  formalism  of  the  singu- 
larity expansion  method.  Well-known  theory  relating  the  asymptotic 
behavior  of  a function  in  the  transform  domain  and  the  support  of  its 
inverse  in  the  time  domain  is  brought  to  bear  on  the  SEM  formulation. 
The  significant  conclusions  are  summarized  in  the  following  list. 

1.  The  support  of  the  time  domain  inverse  kernel  to  an  electro- 
magnetic integral  equation  is  dependent  on  the  asymptotic 
behavior  of  the  transform  domain  form  of  the  inverse  kernel. 
This  asymptotic  behavior,  and  hence  the  support,  depend,  in 
general,  on  the  source  and  observation  points. 

2.  The  dependence  of  the  support  on  the  source  and  observation 
points  affects  the  coupling  coefficient  form  applicable  to 
the  formulation  for  a particular  problem. 

3.  In  a numerically  treated  SEM  problem,  we  have  only  a limited 
amount  of  pole  data  available.  These  data  are  adequate  to 
form  the  components  of  an  approximate  solution  by  way  of  the 
transform  domain  inverse  kernel.  It  is  not  adequate  to  pro- 
vide asymptotic  behavior  of  the  inverse  kernel. 
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4.  We  must  turn  to  other  information  in  order  to  gain  this  asymp- 
totic information  which  is  necessary  to  construct  correctly  the 
time  domain  SEM  solution,  in  particular,  the  coupling  coefficients. 

5.  The  observation  of  a numerically  computed  time  domain  inverse 
kernel  allows  us  to  choose  a coupling  coefficient  form  consistent 
with  causality. 

By  example,  it  is  seen  that  for  at  least  one  problem,  that  of  a straight- 
wire  scatterer,  direct  application  of  the  class  1 coupling  coefficient  form 
introduces  acausal  contributions  into  the  solution.  This  observation  is 
consistent  with  the  experiences  of  Tesche  and  of  Wilton  and  Umashanker  in 
attempts  to  apply  class  1 coefficients  to  numerical  SEM  solutions  [12],  [13]. 
Further,  Van  Blaricum  and,  independently,  Poggio  have  observed  that  in 
applying  Prony's  method  to  electromagnetics  an  exponential  series  alone 
cannot  describe  the  time  dependence  of  observed  waveforms  while  the  driving 
field  is  present  [14],  [15].  This  observation,  in  effect,  recognizes  a 
breakdown  in  the  class  1 coupling  coefficient  SEM  form  for  these  problems. 

The  class  1 form  relegates  all  time  dependence  to  the  exponential  factors 
in  the  singularity  expansion.  It  indicates  that  other  terms  in  the 
representation  must  indeed  be  time  dependent.  This  too,  is  consistent 
with  the  observation  in  this  chapter  that,  at  least  for  some  problems,  a 
coupling  coefficient  which  varies  with  time  during  the  passage  of  the 
incident  wave  is  required.  A second  potential  source  of  the  additional 
time  variation  is  an  entire  function  contribution  during  the  passage  of 
the  wave. 

The  discussion  in  this  chapter  is  based  on  the  objective  of  writing 
the  SEM  representation  for  a time  domain  solution  in  terms  of  the  residue 
series  alone.  We  observe,  by  example,  that  for  some  structures,  the  class  1 
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coupling  coefficient  introduces  an  acausal  contribution  to  the  solution. 

Baum  [10]  observes  correctly  that  in  the  Mittag-Lef f ler  representation  for 

■-1 

a meromorphic  function  an  entire  function  with  zeros  at  the  poles  of  F 
may  be  added  and  subtracted  so  as  to  modify  the  asymptotic  behavior  of  the 
pole  series  portion  of  the  representation.  The  class  1 form  results  for 
the  series  by  such  a manipulation,  viz.,  from  (2.9)  and  (2.13)  write 


F 1 (r , r ' , s)  = l 6 (s  - s ) 1 M (r)  C (r')  + r ^(r.r'.s) 
m m m m e 

m 


(2.44) 


where  T ^ is  an  entire  function  in  s.  Then  add  and  subtract  an  entire 
e 

function  with  zeros  at  s = s to  obtain 

m 


=-l v -s(t'-p‘r')  r , v-1 

‘ (r,r  ,s)  = e H \ Bm(s  - sm) 


s (t ’-p* r* ) 


M (r)  C (r') 
m m 


+ “ei (r,r' , s ) 


r£1(r,r' ,s)  + 'e^(r,r' ,s)  , 


(2.45) 


where  f , is  the  modified  entire  function.  Now  F ^(r,r',s)  is  a unique 
function  and  its  asymptotic  form  is  the  same  whether  it  be  represented 
by  (2.44)  or  (2.45).  However,  if  we  invert  the  product  of 

1 (r ,r' ,s) , E(r' ,s))  from  (2.45)  we  have  the  superposition  of  the 
series  contribution  and  the  entire  function  contribution.  Then,  treating 
the  series  inversion  termwise,  we  effectively  have  introduced  the  class  1 
asymptotic  exponent,  T = p • r'  + t',  to  each  term  in  the  series. 

la 

« 

Thus,  a class  1 form  results  for  the  series  in  the  time  domain: 

J(r,t)  - Jj.U.t)  + Jg,  (r ,t)  , 


(2.46) 
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with  a current  expanded  by  the  class  1 expansion  of  (2.36)  and  (2.39). 
The  examples  of  the  transmission  line  demonstrate  that  potentially 
includes  acausal  contributions.  Therefore,  the  manipulation  between 
(2.44)  and  (2.45)  necessarily  introduces  a contribution  in  Je,(r,t) 
which  cancels  the  acausal  portion  of  J„.  The  following  one-term  example 
illustrates  the  principle.  Let 


f(s) 


f(t)  = u ( t ) e 


V 


The  partitioning  below  results  in  a sum  form: 

/(s  - Sq) 

The  second  term  is  an  entire  function.  Termwise  inversion  of  this  ex- 
pression yields 

g £ g £ St 

f(t)  = u(t  + Tq)  e ° - e ° [u(t  + TQ)  - u(t) ] = u(t)  e ° 

In  this  example  f(s)  is  the  transform  of  a "causal"  function.  Modifying 
the  asymptotic  form  pole  term  allows  the  premature  application  of  its 
inverse  transform  but  at  the  expense  of  an  additive  canceling  term. 

By  considering  the  net  asymptotic  form  of  V ^(r,r',s)  as  in  the 
preceding  sections,  we  obtain  a time  domain  representation  involving 
only  a pole  series  but  with  a coupling  coefficient  of  the  class  2 or 
class  X type  which  is  time  dependent.  By  partitioning  the  asymptotic 
behavior  between  the  series  and  entire  functions  as  in  (2.45),  we  gain 
a class  1 pole  series  form  at  the  expense  of  a required  entire  function 
to  cancel  any  acausal  behavior  present  in  the  class  1 pole  series. 


f (s) 


-Vvs) 


/(s  - 


s0}  + 


1 - e 


"To(so"s) 
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3.  DETERMINATION  OF  THE  SEM  REPRESENTATION 
FOR  A RECTANGULAR  APERTURE 


3 . 1 Introduction 

This  chapter  describes  the  numerical  determination  of  the  SEM  quantities 
of  poles,  natural  modes,  coupling  vectors  and  normalization  constants. 

The  determination  is  based  on  a method  of  moments  numerical  solution  to 
the  integrated  E-field  integral  equation  form  for  the  dual  problem  of  a 
rectangular  plate  [1].  These  coupled  integral  equations  have  a form 
counterpart  to  the  Hallen  integral  equation  for  a linear  antenna. 

In  Section  3.2  we  develop  the  integral  equations  for  the  problem 
for  complex  frequencies  and  develop  the  symmetry  relationships  which 
exist  between  the  natural  modes  for  the  problem.  The  symmetry  provides 
a significant  computational  benefit  in  the  numerical  solution  of  the 
integral  equation. 

Section  3.3  describes  tlv  numerical  solution  procedure  based  on  the 
method  of  moments  used  in  solving  the  integral  equations.  Other  more 
routine  computational  considerations  are  discussed  there,  too.  In 
attacking  a two-dimensional  SEM  problem  numerically,  it  is  essential 
to  be  mindful  of  computational  efficiency. 

Section  3.4  presents  the  results  of  some  numerical  checks  used  to 
validate  the  solution  procedure  and  the  programming  thereof.  These 
checks  include  convergence  studies  for  the  poles. 

Section  3.5  describes  the  pole  search  strategy  used  in  obtaining 
the  SEM  data  for  the  aperture.  The  use  of  argument  number  in  pole  lo- 
cation is  discussed  from  an  applications  point  of  view. 
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3.2  Thln-Plate  Integral  Equation  Formulation  for  Complex  Wavenumber 
3.2.1  The  integral  equations  for  complex  frequency 
Rahmat-Samli  and  Mittra  [1]  give  an  integral  equation  formulation 
for  the  rectangular  plate  subject  to  time-harmonic  excitation.  Their 
results  may  be  directly  extended  to  the  complex  wavenumber  case.  That 
is,  for  the  geometry  in  Figure  3.1  with  exp[st]  time  dependence, 
s = a + jo)  complex,  and  an  incident  plane-wave  magnetic  field  component 


H1  = [HoXax  + Hoyay  + HoZV  exp[sp  • r/c],  the  following  coupled 
integral  equations  result: 


y2  y2(jx(x’ y)i 

-L/2  -w/2  |J  (x,y) j 


r H 

K(x,y  | x'  ,y ' ) dx'  dy'  = / °^exp[px  x + py  y)/c] 

Z ljHox> 


J I Cntj1^1  exp [ j (n  + 1)$]  J^^-jsp/c) 


+ 


1 


exp [ j (n 


1)4>]  Jn_1(-jsp/c)  ] 


(3.1) 


The  kernel  is  given  by 

K(x,y|x',y')  - exp[-sR/c]/R 


(3.2) 


with 


R » l(x  - x')2  + (y  - y’)2]1/2 


The  propagation  is  in  the  direction  of  the  unit  vector  p.  The  J (x,y) 


and  J^(x,y)  denote  the  respective  x and  y components  of  current  on  the 

plate;  J (s)  denotes  the  Bessel  function  of  the  first  kind;  C are  un- 
n n 

known  constants;  c is  the  velocity  of  light;  and  (p,4>)  are  the  polar 
coordinates  for  the  point  (x,y)  on  the  plate.  Equation  (3.1)  holds 
for  x e (-L/2,L/2)  and  y e (-w/2,w/2),  and  z = 0. 

It  is  pointed  out  that  the  tvs  integral  equations  represented  by 
(3.1)  are  coupled  through  the  in  the  summation  in  the  right-hand  side 
This  summation  is  simply  a Bessel  function  expansion  of  the  homogeneous 
solution  to  the  wave  equation  which  occurs  in  the  derivation  of  (3.1). 

The  current  solutions  to  (3.1)  satisfy  the  Meixner's  edge  condition 
[2];  namely. 


J 

x 


J 

y 


[±  (L/2  - d) , y]  -+•  d 
[±  (L/2  - d)  , y ] -*■  d 


1/2  ^ 


Jx[x,±(w/2  - d)] 
Jy[x,±(w/2  - d)] 


- d"1'2 

- d1/2 


0 


(3 


The  first  and  fourth  of  these  require  that  the  normal  component  of  cur- 
rent vanish  at  the  edges  of  the  plate.  It  is  convenient  to  enforce 
these  conditions  directly  in  the  numerical  procedure  as  a means  of  eval- 
uating the  constants  C^.  This  procedure  is  described  more  fully  in 
Section  3.3. 


3.2.2  Symmetry  conditions  for  the  natural  mode  currents 

The  natural  frequencies  of  (3.1)  occur  when  the  complex  frequency 

s is  such  that  there  are  nontrivial  J and  J and  the  accompanying  C 

x y n 
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which  satisfy  (3.1)  for  H =0.  Such  J and  J solutions  are  natual  mode 

x y 

current  solutions  for  the  rectangular  plate,  and  the  concomitant  value 
of  s is  a pole  of  the  plate.  The  vanishing  of  incident  wave  dependence 
gives  rise  to  symmetry  in  the  integral  equations.  By  discerning  the 
symmetry  relations  a priori  and  bringing  them  to  bear  upon  solution  pro- 
cedures, one  gains  significant  computational  savings  in  the  numerical  solu- 
tion for  poles  and  natural  modes. 

The  excitation-free  form  of  (3.1)  is 


L/2  w/2  j cc  , 

/ / J K(x,y|x',y’)  dx’  dy'  = ^ \ C (jn+1 

-L/2  -w/2  X s _co  n\ 


exp [ j (n  + 1)4>]  Jf^fl(-jsp/c) 


+jn  1 exp [ j (n  - 1)<J>]  Jn_1(-jsp/c)|  (3.4a) 


L/2  w/2  co  s' 

/ / Jv  K(x,y|x',y')  dx'  dy ' = £ C <jn  exp[j(n  + 1)4>]  J (-jsp/c) 

-L/2  -w/2  y s -cc  n n+i 


j exp [ j (n  - 1) 4> ] (-jsp/c)>.(3.4b) 


By  using  the  symmetry  of  the  Bessel  function  with  respect  to  order, 
expanding  the  exponentials  by  way  of  Euler's  identity,  and  appropriately 
adjusting  the  indices,  one  arrives  at  the  following  equations  after 
some  manipulation. 


L/2  w/2  oo  r 

I / Jx  K dx*  dy'  = 4^  1 )Jn+  dn  [ cos (n  + 1)<J>  J (-jsp/c) 

-L/2  -w/2  x s n*0l  " n+1 


u . cos  (n  - 1)<1>  J . (-jsp/c)]  - jn  d [sin  (n  + 1)<J>  J (-jsp/c) 
ii“  i n- 1 n tvri 


un_1  sin(n  - 1)$  Jn_1  (-jsp/c)j 


(3.5a) 
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and 


L/2  w/2 

J / J K dx'  dv' 
-L/2  -w/2  y 


+ u , sin(n  - l)<i)  J 
n-I  n- 


7TC 

S 


l 

n=0 


1(-jsp/c)  ] 


jn+1  d+[sin(n  + 1 ) <J>  Jn+1(-jsp/c) 

+ jn  dn [ cos  (n  + 1) <t>  Jn+1(-jsp/c) 

+ un-1  cos  (n  - 1) 4>  J ^ (-jsp/c)  ] 


, (3 


where 


d"  = C ± C 
n n -n 


and 


I1’ 

n > 0 

l°- 

n < 0 

It  is  noted  that  the  d^  multiply  terms  containing  cosine  functions  in 

the  equation,  while  they  multiply  terms  containing  sine  functions 

in  the  J equation.  The  situation  is  reversed  for  the  d . 

y n 

Because  of  the  symmetry  properties  of  the  kernel,  the  integral 
operator  on  the  left-hand  sides  of  (3.5)  produces  a function  whose 

symmetry  character  is  identical  to  that  of  the  current  on  which  it 

+ 

operates.  Then,  for  a given  current  symmetry,  only  part  of  the  d^  on 
the  right-hand  side  may  be  nonzero  because  of  the  symmetries  possessed 
by  the  trigonometric  terms.  Thus,  the  respective  symmetries  for 
and  Jy,  which  are  compatible,  and  the  surviving  terms  in  the  right-side 
series  may  be  discerned  by  (1)  postulating  a symmetry  for  J , (2)  deter- 
mining from  (3.5a)  which  right-hand-side  terms  survive  so  as  to  be 


.5b) 
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compatible  with  the  symmetry,  (3)  observing  in  (3.5b)  the  variation 
of  the  terms  in  the  right-hand  side  which  have  nonzero  coefficients, 
and  (4)  determining  the  symmetry  conditions  from  the  symmetry  ob- 
served in  the  right-hand  side  and  thus  compatible  with  the  postulated 
J symmetry  conditions. 

For  example,  if  is  symmetric  with  respect  to  the  y-axis  and 
antisymmetric  with  respect  to  the  x-axis,  only  sin(n  + l)cj>  terms  with 
n even  are  compatible  in  (3.5a).  Thus,  only  d^,  n even,  may  be  nonzero. 

In  the  right-hand  side  of  (3.5b),  the  coefficients  multiply  cos (n  + 1)$ 
terms  with  n even.  These  cosines  sum  to  functions  which  are  antisymmetric 
with  respect  to  the  y-axis  and  symmetric  with  respect  to  the  x-axis. 

Stated  mathematically,  if 


Jx(x,y)  = Jx(-x,y) 

and 

Vx.y>  = -Jx(x>~y> 

then 

d+  - 0, 
n 

d“  = 0, 
n 

and 

JyU.y)  - -J(-x,y) 
Jy(x»y)  - Jy(x.-y) 


for  all  n , 
n odd  , 


(3.6a) 


(3.6b) 


(3.6c) 

(3.6d) 


(3. 6e) 
(3. 6f ) 
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These  vector  symmetries  are  in  accord  with  the  general  symmetry  relations 
given  by  Baum  [3].  The  information  in  (3.6)  may  be  used  to  reduce  the 
complexity  of  the  integral  equations  (3.4),  viz.,  by  (3.6a,b,e,f)  the 
range  of  each  integration  may  be  halved  while  by  (3.6c,d)  the  zero  terms 
of  the  right-hand  side  are  known  a priori: 

L/2  w/2  » 

I / J K (x,y|x',y')  dx'  dy'  = — ^ £ d jn  i[sin(n  + 1)$  J (-jsp/c) 

0 0 X S n=0  n n_i 

n even 

- sin (n  - 1)4)  J .(-jsp/c)]  (3.7a) 

n-i 


and 


L/2  w/2 


/ / J K+  (x,y | x' ,y ' ) dx ' dy'  = ~ £ jn+i  d [cos(n  + 1)$  J (-jsp/c) 


. _ JL9l  T ,n+l  .,-1 


0 0 


n-0 
n even 


n-»-l 


+ un-1  cos  (n  - 1) <f>  Jn_1(-jsp/c)  ] (3.7b) 


where 

K+  (x,y | x' ,y ' ) = K(x,y|x',y')  - K(x,y | -x' ,y ' ) 

+ K(x,y | x' ,-y ' ) - K(x,y | -x' , -y ' ) (3.8a) 

and 

K_+ (x,y | x' ,y ' ) = K(x,y|x',y')  + K(x,y | -x' ,y ' ) 

- K(x,y|x,,-y')  - K (x,y | -x' ,-y ' ) . (3.8b) 


For  subsequent  reference 


(x,y  | x'  ,y ' ) = K(x,y|x’,y')  + K(x,y ) -x' ,y' ) 

+ K(x,y| x' ,-y’)  + K (x,y j -x' ,-y ' ) 

and 
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(3.8c) 


K (x,y | x’ ,y ' ) = K(x,y|x',y')  - K(x,y |-x' ,y ' ) 

- K(x,y | x' ,-y ' ) + K(x,y|-x' ,-y’)  (3.8d) 

are  defined  as  well.  Equations  (3.7)  are  enforced  for  z = 0,  x e (0,L/2) 
and  y e (0,w/2). 

Table  3.1  summarizes  the  four  symmetry  cases  which  are  derived  as  in 
the  foregoing  discussion.  By  means  of  this  table,  four  integral  equation 
pairs  can  be  constructed  in  the  spirit  of  (3.7)  by  replacing  the  kernels 
in  (3.7^  with  the  appropriate  kernels  from  the  table  and  retaining  only 
t'ie  nonvanishing  terms  in  the  series  expansion. 

Figure  3.2  depicts  qualitatively  the  respective  modal  current  dis- 
tributions for  the  lowest  frequency  natural  resonance  exhibiting  each 
symmetry . 

3.3  The  Numerical  Model 

3.3.1  Discretization  of  the  integral  equations 

The  integral  equation  pair  of  the  form  (3.7)  for  each  of  the  four 
symmetry  cases  can  be  discretized  by  the  method  of  moments.  In  the 
work  reported  here,  two-dimensional,  subsectionally  constant  expansion 
functions  were  used  with  collocation  testing.  The  zoning  scheme  is 
represented  in  Figure  3.3. 

The  unknown  currents  and  were  expanded  in  piecewise  constant 
functions  as  in  [1]  with  subsectioning  of  the  form  given  in  Figure  3.3. 
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Notice  that  half-width  patches  are  used  at  the  edges  of  the  plate  so 
that  match  points  lie  precisely  on  the  edge  of  the  plate.  The  half- 
width pulse  has  proved  useful  in  realizing  the  actual  electrical  size 
of  a body  in  one-dimensional  problems  [4]. 

The  boundary  condition  J “0  must  be  enforced  on  selected  patches 
1 norm  r 

at  the  edge  of  the  plate  as  discussed  in  [1].  Concomitantly,  only  as 

■f  . 

many  d^  s are  retained  in  the  right-hand-side  summation  in  (3.7)  as 
there  are  current  values  preassigned  to  zero.  The  shaded  patches  in 
Figure  3.3  indicate  the  selection  of  patches  where  a current  component 
is  preassignea  a zero  value.  At  the  corner  patch,  both  components  are 
preassigned  zero  values.  By  assigning  one  match  point  per  expansion 
patch  and  by  retaining  one  series  expansion  term  for  each  current  value 
preassigned  in  each  of  the  two  integral  equations,  we  obtain  a square 
system  of  linear  equations.  The  truncated  summation  is  taken  to  the 
left-hand  side  so  that  a homogeneous  system  results.  The  matrix  organ- 
ization used  to  represent  these  equations  is  given  by  the  partitioned 
matrix  equation 

k(8>]  H I 1 K) 

HjW  [-0  - 

[°  ] H L 1 F 


The  symbols  are  as  follows.  The  matrix  [M^]  is  the  moment  matrix  for 
the  integral  equation  as  given  in  (3.7a)  but  with  the  kernel  symmetry 
appropriate  to  the  mode  symmetry  for  which  an  individual  solution  is 
being  conducted.  The  matrix  [My]  is  the  Jy  moment  matrix  as  in  (3.7b) 
with  the  kernel  symmetry  compatible  with  that  of  the  kernel.  The 
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coupling  enters  through  [M^ ] - This  matrix  is  formed  by  the  negative  of 
sine-Bessel  function  products  on  the  right-hand  sides  of  (3.7a,b)  which 
survive  under  the  symmetry  considered.  Sufficient  terms  are  taken  in 
the  series  to  render  a square  system  after  entries  for  the  edge  currents 
which  are  preassigned  to  zero  are  deleted  from  moment  matrices. 

In  using  this  integral  equation  in  the  real  frequency  domain, 
Rahmat-Samii  and  Mittra  experienced  a numerically  unstable  matrix  if  the 
ro  current  edge  condition  was  enforced  at  every  zone  on  the  edge 
of  the  plate.  This  difficulty  was  attributed  to  redundancy  in  the 
imposition  of  the  boundary  condition  and  was  corrected  when  the  con- 
dition was  enforced  only  for  every  second  or  third  zone  at  the  edge. 

The  edge  current  zones  in  between  which  were  left  undetermined  attained 
a value  of  zero  within  limits  of  numerical  error  in  the  course  of  the 
numerical  solution. 

The  procedure  of  preassigning  currents  at  preselected  edge  zones 
proved  to  fail  when  it  was  applied  in  the  SEM  pole  search  in  the  complex 
s-plane.  It  was  observed  that  the  moment  matrix  had  zeros  in  the  com- 
plex plane  whose  associated  homogeneous  solution  did  not  approximate 
zero  for  the  zones  in  between  specified  patches.  The  procedure  of  spec- 
ifying at  every  edge  zone  that  the  normal  current  component  vanish 
provided  physically  consistent  solutions. 

The  matrix  in  (3.9)  is  a function  of  the  complex  frequency  s. 

An  SEM  "pole"  or  complex  resonance  occurs  when  this  composite  matrix 
has  a zero  determinant.  It  follows  directly  from  the  analyticity  of 
the  kernels  of  the  original  integral  equations  (3.2)  that  the  determinant 
is  an  analytic  function  of  s throughout  the  whole  of  the  complex  s- 
plane.  The  Muller  method  for  searching  for  roots  of  nonlinear  systems 
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of  equations  provides  an  efficient  algorithm  [5]  for  systematically 
searching  for  the  values  of  s for  which  the  determinant  vanishes. 

Many  evaluations  of  the  determinant  are  required  in  the  course  of  the 
iteration  to  locate  a zero.  This  places  a premium  on  computational 
efficiency  in  the  evaluation  of  the  matrix  and  its  determinant.  The 
next  subsection  discusses  the  means  taken  to  attain  computational  ef- 
ficiency. 


3.3.2  Algorithmic  considerations  in  evaluating  the  system  determinant 

Some  considerations  taken  into  account  in  generating  the  system 
matrix  and  evaluating  its  determinant  efficiently  are  discussed  in  this 
section.  Since  these  two  operations  must  be  repeatedly  carried  out  for 
many  values  of  s in  the  course  of  determining  the  natural  frequencies  of 
the  plate,  it  is  essential  that  clean,  efficient  computer  programming 
and  coding  be  used  so  that  execution  of  the  program  will  be  affordable. 

The  volume  of  code  in  the  algorithms  is  consistently  compromised  toward 
a larger  size  in  order  to  meet  the  following  two  time-efficient  objectives: 

1.  Avoidance  of  calculating  the  same  quantity  twice;  and 

2.  Avoidance  of  logical  decisions,  particularly  those  which  might 
be  imbedded  in  loops. 

The  program  is  discussed  in  the  context  of  the  following  major 
segments : 

1.  Computation  of  an  "interaction  matrix"; 

2.  Construction  of  the  nonzero  submatrices  of  the  system  matrix 
from  the  interaction  matrix; 

3.  Computation  of  the  series  terms'  submatrix;  and 


4. 


Determinant  evaluation. 
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The  major  contribution  to  the  elimination  of  redundant  calculations 
is  the  one-time  computation  of  an  "interaction  matrix"  which  is  made  up 
of  the  individual  kernel  integral  terms  from  (3.2)  for  all  argument 
combinations  which  occur  in  the  computation.  The  subsequent  program 
step  then  picks,  by  subscript,  entries  from  this  matrix  and  constructs 
the  appropriate  kernel  from  one  of  Equations  (3.8)  according  to  the 
symmetry  conditions  being  solved.  This  procedure  can  be  viewed  in  terms 
of  the  layout  given  in  Figure  3.4a.  The  terms  in  the  interaction  matrix 
are  those  evaluated  for  the  match  point  as  shown  in  the  lower  left 
with  the  source  patches  indexed  over  the  entire  plate  to  generate  the 
matrix.  Thus,  all  geometric  relationships  which  occur  in  the  kernel 
terms  are  encompassed  in  the  calculation.  Note  that  all  source  patches 
are  full  patches  for  this  calculation.  The  effect  of  half  patches  at 
the  edges  is  accounted  for  by  weighting  by  a factor  of  1/2  the  edge 
contributions.  The  kernel  integral  appropriate  to  the  symmetry  is 
constructed  by  summing  with  correct  signs  the  appropriate  elements  from 
the  matrix.  Figure  3.4b  gives  an  example  of  the  four  source  patches 
entering  into  one  kernel  integral. 

Differing  degrees  of  sophistication  are  required  in  the  calculation 
of  the  interaction  terms  depending  on  the  spacing  of  the  patches  for 
which  an  interaction  is  being  calculated.  For  the  self  patch,  i.e., 
the  patch  in  which  the  match  point  resides,  the  integration  of  the  kernel 
must  be  performed  analytically  because  of  the  integrable  singularity  in 
the  kernel  there.  For  the  patches  adjacent  to  the  patch  containing 
the  match  point,  the  kernel  is  a rapidly  varying  but  well-behaved  func- 
tion. The  integration  over  these  patches  is  evaluated  numerically  by 
a polynomial  approximation.  For  patches  further  separated,  the  kernel 
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is  slowly  varying  and  the  integral  is  evaluated  approximately  as  the 
product  of  the  value  of  the  kernel  at  the  center  of  the  patch  and  the 
area  of  the  patch. 

Some  minor  time  economy  is  achieved  in  the  matrix  filling  algorithm, 
which  is  a four-dimensional  loop.  The  economy  is  found  in  the  form  of 
decision-free  indexing,  that  is,  the  source  contributions  from  interior 
patches,  from  |x|  = L/2  edge  patches,  from  |y|  = w/2  edge  patches,  and 
from  corners  take  on  different  forms.  Rather  than  index  over  all  source, 
patches  with  logical  decisions  implemented  to  discriminate  among  the 
four  cases  above,  four  different  loops  are  used. 

The  computation  of  the  series  term  submatrix  is  relatively  straight- 
forward. Because  the  Bessel-trigonometric  products  appear  in  two  terms 
each,  they  are  all  precalculated  and  stored  in  a vector.  The  individual 
terms  are  then  constructed  from  them. 

The  determinant  evaluation  profits  significantly  from  an  exploitation 
of  the  sparseness  of  the  matrix.  Either  of  two  approaches  may  be  taken 
to  the  sparse  matrix  manipulations.  One  is  to  separate  the  matrix 
algebraically  and  calculate  an  inverse  as  a composite  of  inverses  of 
terms  involving  the  submatrices.  The  alternative  approach  is  to  attack 
the  matrix  directly  with  Gaussian  elimination  exploiting  the  sparse- 
ness directly  in  the  algorithm.  The  latter  approach  was  chosen  for  the 
present  purpose  because  it  is  judged  to  be  slightly  faster  computationally 
and  because  in  order  to  determine  natural  mode  solutions  for  the  SEM 
formulation,  the  homogeneous  system  of  equations  occurring  at  a pole 


must  be  backsolved. 
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3.3.3  Numerical  evaluations  of  natural  modes,  coupling  vectors, 
and  normalization  constants 

In  addition  to  the  poles,  the  associated  coupling  vectors  and 
natural  mode  vectors  defined  by  (2.11)  and  the  normalization  constants 
defined  by  (2.13)  form  the  SEM  description  of  a given  problem.  In  a 
method  of  moments  solution,  these  are  calculated  in  a discretized 
fashion  and  are  "one-time"  calculations  for  each  pole.  We  need  not  be 
quite  so  mindful  of  computation  costs  as  we  are  in  the  matrix/determinant 
computation  which  is  used  iteratively. 

The  natural  mode  computation  is  a straightforward  one  once  a pole 
is  located.  The  Gaussian  elimination  with  maximum  pivoting  used  to 
calculate  the  determinant  in  the  course  of  pole  location  results  in  a 
triangularized  matrix  which  has  a zero  in  the  lower-right  position  when 
a zero  determinant  is  computed.  Thus,  the  last  triangularized  matrix 
upon  return  from  the  zero  search  iteration  may  be  backsolved  by  assigning 
a value  to  the  last  element  in  the  vector  of  unknowns  and  then  backsolving 
the  system  of  equations  systematically.  The  present  work  applies  this 
procedure  using  the  sparse  matrix  format  so  as  to  be  compatible  with 
the  determinant  evaluation  algorithm  discussed  above. 

The  coupling  vectors  are  found  in  a similar  manner  except  that  the 
transpose  conjugate  of  the  moment  matrix  must  be  triangularized  and 
backsolved  for  the  adjoint  solution.  The  transpose  conjugate  is  that 
of  Equation  (3.9).  The  benefit  of  sparseness  is  not  as  great  for  the 
transposed  form  because  the  coupling  there  appears  in  the  lower  rows 
of  the  matrix.  The  moment  matrix  can  be  recalled  from  storage  or  re- 
computed and  the  transpose  conjugate  formed  in  the  composite  format 
of  (3.9).  This  matrix  is  triangularized  and  backsolved  with  elementary 
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3 . A • 2 Pole  convergence  for  larger  aspect  ratios 

The  thln-strip  tests  on  convergence,  of  course,  tell  us  very  little 
about  the  formulation  for  apertures  of  larger  aspect  ratios.  The  trans- 
verse component  of  current  is  quite  small  on  the  strip,  and  the  trans- 
verse variation  of  the  longitudinal  component  is  not  significant  over 
the  small  electrical  width  of  the  strip. 

Figure  3.6  provides  some  insight  into  the  convergence  of  the  poles 
for  a range  of  aspect  ratios  for  several  zoning  densities.  The  pole 
trajectory  is  for  the  second  ++  (w.r.t.  J^)  in  the  "layer"  nearest  the 
imaginary  axis.  We  note  that  wL/cir  = 3.0  corresponds  to  a 1-1/2 
wavelength  structure  so  that  a density  of  five  zones  per  L/2  pulses  in 
the  first  quadrant  is  equivalent  to  6.67  pulses  per  wavelength. 

The  convergence  of  pole  locations  for  larger  aspect  ratios  was 
tested  by  comparing  those  obtained  with  4 x A zones  in  the  first  quadrant 
with  those  computed  with  5x5  zones . The  comparisons  are  shown  as 
the  solid  and  dashed  lines .respectively , in  the  "loop"  region.  The 
aspect  ratios  are  shown  as  decimal  fractions.  We  observe  excellent 
agreement  in  the  two  trajectories.  The  error  at  w/L  = 1.0  is  1.16  per- 
cent of  the  magnitude  of  the  poles. 

In  following  a pole  trajectory,  it  is  desirable  to  use  the  same 
number  of  zones  insofar  as  possible  so  that  a smooth  trajectory  is 
determined.  This  is  not  always  possible  for  smaller  aspect  ratios 
because  of  numerical  instability  problems  inherent  in  over-dense  zoning. 
The  continuation  of  the  A x A zone  trajectory  is  shown  down  to  w/L  = 0.3 
in  Figure  3.6.  For  aspect  ratios  smaller  than  0.3,  numerical  instability 
besets  us  and  the  trajectory  cannot  be  determined.  In  practice,  the 
numerical  instability  may  be  manifested  in  one  of  two  forms:  either 
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a failure  to  converge  in  the  pole  search  or  an  abrupt  "tailing  off"  in  a 
smooth  trajectory.  The  former  manifestation  forced  the  termination  of  the 
4x4  trajectory.  The  latter  is  evident  in  the  4x3  trajectory  plotted. 
The  filamentary  dipole  pole,  which  is  at  sL/cir  = j 3.0,  is  the  limiting 
point  for  this  trajectory,  and  the  4x3  result  is  clearly  departing  from 
this.  The  4x2  trajectory  is  seen  to  be  consistent  with  the  correct 
behavior.  Therefore,  in  steps  of  0.1  in  aspect  ratio,  we  might  choose 
the  4x4  trajectory  between  1.0  and  0.3,  the  4 x ^ point  for  0.2,  and 
the  4x2  point  for  0.1  to  obtain  a smooth,  physically  consistent  tra- 
jectory. 

This  strategy  was  followed  in  obtaining  the  pole  trajectories 
given  in  the  next  chapter.  Namely,  use  4x4  zones  beginning  with  the 
larger  aspect  ratios  on  a trajectory;  follow  the  trajectory  as  the 
aspect  ratio  decreases,  reducing  the  transverse  zoning  density  only  as 
much  as  necessary  to  obtain  smooth  physically  consistent  results. 

Another  divergence  phenomenon  is  encountered  in  tracing  the  pole 
trajectories  for  -+  and  — modes.  We  observe  that  these  modes  are 
antisymmetric  in  the  narrow  direction  on  a thin  strip.  Thus,  these 
symmetry  conditions  impose  180  degree  phase  opposition  upon  closely 
spaced  currents  on  the  thin  strip.  We  might  anticipate  an  instability 
in  the  integral  equations  themselves  for  these  modes  in  thin  strips. 

(This  instability  is  the  counterpart  to  that  which  one  encounters  when 
solving  two  parallel  linear  antennas  which  are  closely  spaced.) 

Instability  does  appear  in  these  solutions.  Its  manifestation  is 
that  for  thin  strips  with  the  longitudinal  component  of  current 
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antisymmetric  in  the  transverse  direction  the  pole  search  converges, 
but  the  point  to  which  the  algorithm  converges  is  weakly  dependent 
(within  a few  percent)  on  its  starting  point.  In  addition,  the  pole 
values  obtained  depend  strongly  on  the  number  of  zones  used.  This  zone 
dependence  is  evident  in  Figure  3.7.  The  result  is  that  none  of  the 
-+  or  — trajectories  are  determinable  by  the  present  means  for  aspect 
ratios  of  below  0.2  or  0.3.  The  dominant  current  component  is  anti- 
symmetric in  the  narrow  direction,  however.  This  results  in  weak  coupling 
of  these  modes  for  smaller  aspect  ratios  so  that  they  are  generally 
negligible . 

3.5  Pole  Location  Strategies 

The  basic  pole  location  method  used  in  this  work  was  the  Muller 
iteration  to  find  zeros  of  the  determinant  of  the  moment  matrix  as  a 
function  of  the  complex  frequency  s as  was  mentioned  in  Section  3.4. 

It  appears  to  be  a common  trait  to  numerical  SEM  formulations  that 
except  for  the  poles  near  the  imaginary  axis  in  the  complex  plane,  the 
contours  of  the  determinant  indicate  the  presence  of  the  zero  determinant 
in  a highly  local  fashion.  The  paragraphs  below  describe  the  methods 
used  in  locating  the  poles  in  an  exhaustive  fashion  without  undue  dif- 
ficulty with  the  local  character  of  the  poles. 

Baum  [6]  has  proposed  using  the  argument  number  of  the  determinant 
of  the  moment  matrix  as  an  indicator  of  the  presence  of  a zero  of  the 
determinant  within  a prescribed  region  of  the  complex  plane.  The  method 
exploits  the  analyticlty  of  the  determinant  as  a function  of  s.  It 
follows  that  the  number  of  zeros  of  the  determinant  D(s)  within  a region 
in  the  s-plaue  bounded  by  a closed  contour  C is  given  by 
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Nc  * V(,{D  (s)  } - — j-  / arg[D(s)]  ds  , 

c 

where  denotes  the  "variation  of  the  argument"  of  D(s)  around  the 
contour  C.  (See,  for  example  [1].)  This  approach  was  used  conveniently 
and  successfully  in  the  present  work  to  localize  a pole  to  a small 

region  before  using  the  iterative  method  to  locate  it  specifically. 

» 

Figure  3.8  illustrates  the  method  used.  The  region  of  interest 
in  the  complex  plane  was  subdivided  into  moderately  sized  quadrants. 

The  argument  number  was  evaluated  around  each  quadrant  from  twenty-one 
values  of  D(s)  along  each  side  of  the  quadrant.  The  contributions  to 
Nj,  from  each  leg  of  each  quadrant  were  evaluated  separately  so  that  only 
alternate  quadrants  needed  to  be  computed.  The  results  of  the  pro- 
cedure for  the  -H-  modes  are  shown  in  Figure  3.8  for  aspect  ratios  of  1.0 
and  0.5.  The  pole  trajectories  eventually  determined  are  included  for 
reference  purposes.  It  is  seen  that  the  method  correctly  localizes  all 
of  the  poles  and  provides  confidence  that  the  search  is  exhaustive. 

Each  quadrant  was  subsequently  quadrasected . By  evaluating  the  argument 
number  on  the  four  sides  of  two  of  the  four  subquadrants,  the  pole  is 
still  better  isolated.  The  Muller  iteration  was  then  begun  at  the 
center  of  the  subquadrant  to  which  the  pole  was  isolated  and  in  every 
case  converged  from  that  point. 

Once  the  two  points  on  a trajectory  corresponding  to  0.5  and  1.0 
aspect  ratios  were  located,  the  trajectory  was  followed  by  changing  the 
aspect  ratio  slightly  and  by  using  the  pole  location  for  the  original 
aspect  ratio  as  a starting  point.  This  method  of  incrementing  w/L 
fails  in  some  cases  when  the  pole  location  moves  rapidly  with  respect 
to  w/L.  The  following  projection  scheme  proved  useful  in  these  cases. 


We  observe  from  Figure  3.2  that  for  w/L  = 1 a ninety-degree  spatial 
rotation  of  the  ++  mode  results  in  the  — mode  configuration  in  the 
original  coordinate  system.  Similarly,  the  +-  mode  rotates  into  itself, 
and  the  — t-  mode  does,  as  well.  This  means  that  a ++  distribution  for 
w/L  = 2.0  is  identical  to  a — distribution  with  w/L  = 0.5.  Similarly, 
the  pole  for  each  case  will  be  identical  when  normalized  to  the  long 
dimension. 

Figure  3.9  shows  the  ++  pole  trajectory  continued  on  past  w/L  = 1.0 
to  w/L  = 8.0.  Note  that  the  scales  for  the  plot  are  normalized  to  L. 

To  renormalize,  for  example,  the  w/L  = 2.0  -H-  pole  to  w,  we  need  to 
multiply  the  pole  value  sL/cir  by  w/L.  Because  of  the  rotation  property, 
the  renormalized  w/L  = 2.0  ++  pole  is  identical  with  the  — w/L  = 0.5 
pole.  This  relationship  is  denoted  graphically  in  Figure  3.9  by  the  line 
segment  passing  through  the  w/L  = 2.0  point  on  the  ++  trajectory  and 
striking  the  w/L  = 0.5  point  on  the — trajectory.  The  length  of  the 
segment  0B  on  this  line  is  twice  0A.  Three  other  specimen  radii  in- 
dicate this  relationship  at  other  points. 

It  is  observed  that  the  — curve  is  rapidly  varying  with  respect 
to  aspect  ratio.  We  can  follow  it  with  ease,  however,  by  computing  the 
++  curve  on  past  w/L  = 1.0  in  the  incremental  fashion  discussed  above 
and  projecting  it  outward  to  form  the  — trajectory.  This  procedure 
was  employed  to  advantage  in  obtaining  the  results  presented  in  the 
next  chapter,  both  by  projecting  ++  onto  — trajectories  and  by  pro- 
jecting +-  or  -+•  trajectories  onto  themselves.  A posteriori  tests  of 
the  projection  corroborated  that  the  correct  pole  had  been  projected. 
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4.  SEM  CHARACTERIZATION  FOR  THE  RECTANGULAR  PLATE 
4. 1 Introduction 

This  chapter  provides  the  pole  and  natural  mode  data  for  the 
dominant  terms  in  the  SEM  characterization  of  the  rectangular  plate. 

By  duality,  the  representation  applied  to  the  rectangular  aperture  as 
well.  The  data  span  the  three  dominant  resonances  for  the  structure  - 
that  is,  a range  of  frequencies  extending,  nominally,  to  where  the  greater 
dimension  of  the  rectangle  is  one-and-one-half  wavelengths.  The  aspect 
ratio  parameter  w/L  is  varied  between  0.1  and  1.0.  Rahmat-Samii  and 
Mittra  have  observed  in  previously  reported  work  [1]  that  for  frequencies 
corresponding  to  roughly  one  wavelength  and  beyond,  the  aperture  becomes 
quite  broadband  so  that  transient  waves  with  higher  frequency  spectral 
content  are  passed  essentially  unmodified.  This  observation  is  supported 
by  the  present  data  which  show  poles  in  proximity  to  one  another  appear- 
ing at  higher  frequencies. 

Representative  natural  modes  for  selected  poles  are  presented 
herein.  Exhaustive  graphical  presentation  of  the  corresponding  natural 
modes  here  is  prohibitively  cumbersome  both  because  of  their  number 
(approaching  two  hundred)  and  because  of  the  futility  of  trying  to 
represent  complex  vector  functions  of  two  variables.  More  complete 
data  are  available  in  a digital  computer  data  base  format,  however  [2]. 

The  data  base  storage  of  SEM  data  for  complex  structures,  such  as  the 
aperture,  is  more  in  accord  with  ultimate  user  needs,  because  a com- 
puter is  almost  certainly  required  in  any  useful  expansion  of  the  data. 

The  coupling  vectors  and  normalization  constants  are  a part  of 
the  data  base  but  are  not  reported  here.  The  qualitative  aspects  of 
coupling  are  discernible  from  the  natural  mode  data. 
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Section  4.2  deals  with  the  pole  locations  for  the  structure  and 
discusses  the  resonance  behavior  of  the  plate  or  aperture  for  the  four 
modal  symmetries.  Section  4.2  presents  some  of  the  associated  natural 
modes  which  have  been  selected  to  be  representative  of  the  mode  structure 
in  general. 

4 . 2 Pole  Trajectories  of  the  Plate/Aperture  as  a Function  of  Aspect  Ratio 

The  locations  of  poles  for  the  rectangular  structure  are  given  in 
Figure  4.1  for  a region  of  the  third  quadrant  of  the  complex  plane  near- 
est the  origin.  The  third  quadrant  data  suffice  because  all  poles  have 
negative  real  parts  and  are  arranged  with  conjugate  symmetry  in  the 
complex  plane.  The  poles  are  normalized  with  respect  to  L/ctt  where  L 
is  the  greater  dimension  of  the  rectangle.  The  small  numbers  adjacent 
to  trajectories  indicate  the  aspect  ratio  w/L  for  the  tick  mark  beside 
which  the  number  stands.  Ticks  are  spaced  at  aspect  ratio  intervals 
of  0.1.  Dashed  lines  indicate  the  trends  of  trajectories  which  are 
not  completely  determined.  Poles  that  are  left  undetermined  on  a known 
trajectory  either  are  judged  negligible  in  contribution  compared  with 
poles  close  to  the  axis  or  in  the  case  of  - and  — modes  are  undeter- 
minable by  the  present  method  for  the  reasons  outlined  in  Section  3.4.2. 

We  reiterate  that  in  the  latter  case  the  coupling  coefficients  will 
be  small  so  that  we  may  neglect  these  poles. 

In  the  following  paragraphs,  we  consider  the  physical  significance 
of  the  pole  behavior  as  shown  in  Figure  4.1.  We  discuss  the  ++  and 
4—  modes  as  a group  and  then  treat  the  -+  and  — modes.  The  former 
group  possesses  the  common  property  of  being  symmetric  about  the  x-axis, 
i.e.,  in  the  transverse  direction.  The  latter  group  is  antisymmetric 


Figure  4.1.  Pole  trajectories  for  the  four  modal  symmetries.  Tick  marks 
are  spaced  at  increments  of  0.1  in  w/L.  Decimal  fractions 
indicate  specific  w/L  ratio  at  selected  ticks. 


in  the  transverse  direction.  Terminology  of  the  form  ++8  is  used: 
this  refers  to  the  upper-left  trajectory  in  the  -H-  mode  pole  plot, 
which  is  labeled  with  the  large  numeral  8.  We  remind  the  reader  that 
the  symmetry  indication  ++  means  that  the  J component  of  current  is 
symmetric  with  respect  to  the  x-axis  and  with  respect  to  the  y-axis . 

The  ++1 , the  +-1,  and  the  ++2  trajectories  may  be  considered  to 
be  the  fundamental  resonances  for  the  structure.  They  are  the  counter- 
part to  the  "layer  1"  poles  for  the  wire  scatterer  as  described  by 
Tesche  [3].  Indeed,  in  the  thin-strip  limit  (w/L  -*•  0),  these  trajectorie 
approach  those  of  the  thin  cylinder.  Both  the  cylinder  and  the  strip 
trajectories  approach  to  the  nonradiating  filamentary  dipole  resonances 
given  by 


sL/cn  = jl, j2,  ... 


The  dominant  poles  are  observed  generally  to  exhibit  decreasing  "Q" 
(quality  factor)  as  the  width  of  the  structure  increases.  An  exception 
to  this  is  the  "loop"  trajectory  for  the  -4-1-2  pole.  This  trajectory 
shows  a peak  in  the  Q at  w/L  = 0.6.  This  phenomenon  is  due  to  the  fact 
that  the  w/L  =0.6  aspect  ratio  admits  to  simultaneous  resonances  in 
the  and  components  of  current.  This  matter  is  discussed  further 
in  the  next  section. 

The  trajectories,  other  than  these  dominant  three,  show  poles  moving 
toward  the  imaginary  axis  with  increasing  aspect  ratio.  That  is,  for 
thin  strips,  the  dominant  poles  alone  dictate  the  principal  rime  d< 
pendence  of  the  scattered  waveform  on  the  structure.  For  wider 
tures,  the  "layer  2"  poles  move  toward  the  imaginary  axis  and  n 

the  scattered  fields  as  well.  For  example,  for  the  square 
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with  w/L  - 1.0,  we  observe  roughly  equal  order  real  parts  for  six  modes. 
These  complex  frequency  data  corroborate  real  frequency  broadbanding 
for  u)L/cir  >_  1 which  Rahmat-Samii  and  Mlttra  observed  in  [2]. 

The  region  of  the  complex  plane  shown  in  Figure  4.1  is  not  large 
enough  to  allow  comparison  with  the  second  layer  poles  for  the  cylinder 
from  [3],  We  note,  in  passing,  that  in  the  few  observations  made, 
there  is  not  an  agreement  between  the  w/L  -►  0.1  limits  of  the  nondominant 
trajectories  in  Figure  4.1  and  the  layer  2 poles  of  the  cylinder.  It 
is  not  at  all  clear  whether  such  an  agreement  should  exist.  If  it 
should,  either  the  failure  to  go  to  the  filamentary  limit  or  simply 
numerical  errors  inherent  in  the  deep-plane  results  preclude  observa- 
tion of  it. 

The  -+  and  — modes  cannot  be  directly  related  to  cylinder  modes. 

This  is  due  to  the  nonphysical  nature  of  antisymmetry  of  the  longitudinal 
current  in  the  transverse  direction  in  the  thin  strip  limit  of  these 
modes.  It  is  difficult  to  separate  out  a dominant  set  of  poles  from 
Figure  4.1  for  these  symmetries.  We  do  observe  the  grouping  of  many 
poles  in  the  square-structure  limit  very  much  in  the  manner  of  the  ++ 
and  +-  modes.  However,  several  trajectories  move  on  a path  roughly 
parallel  to  the  imaginary  axis.  This  increase  in  Q with  the  narrowing 
of  the  structure  is  thus  an  increase  in  oscillation  frequency  of  these 
modes  at  an  approximately  constant  damping  rate. 

We  make  one  other  final  observation  about  the  poles.  It  is  ob- 
served in  Section  3.5  that  on  a square  structure  the  ++  and  — modes  are 
identical  with  a 90  degree  rotation  of  coordinates.  A consequence 
of  this  property  is  that  there  is  a ++  w/L  ■ 1.0  pole  identical  with 
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each  — w/L  • 1.0  pole.  This  is  evident  in  Figure  4.1.  For  example, 
the  — 1 and  ++1  poles  are  Identical  for  w/L  ■ 1.0.  The  numbering 
scheme  for  the  trajectories  is,  in  fact,  chosen  so  that  the  -H-n  and 
— n trajectories  coalesce  at  w/L  - 1.0  for  each  trajectory. 

4 . 3 Natural  Modes  for  the  Rectangular  Structure 

4.3.1  Graphical  representation  of  the  modes 

In  the  four  subsections  which  follow,  we  present  graphical  data 
Intended  to  characterize  the  natural  modes  associated  with  the  poles 
shown  in  Figure  4.1.  Complete  mode  data,  i.e.,  data  for  each  pole 
increment,  are  precluded  by  the  sheer  volume  of  Information  involved. 

As  a result,  representative  mode  data  selected  from  the  complete  data 
base  are  given  in  the  figures  which  follow. 

Each  natural  mode  is  a two-component  complex-valued  vector  function 
of  two  variables.  The  form  of  the  displays  is  evident  in  Figure  4.2. 
Each  mode  is  plotted  as  a current  distribution  on  a thin-plate  scatterer 
and  separate  plots  for  the  x and  y current  components  are  given.  The 
curves  are  magnitude/phase  representations  of  the  complex  quantities. 

The  rectangular  shape  of  the  plate  is  shown  in  proper  proportion  for 
each  of  the  current  components . The  graph  below  the  rectangle  gives  the 
appropriate  current  component  along  a cut  in  the  transverse  coordinate 
y.  The  graph  to  the  right  gives  the  same  current  component  along  a cut 
in  x.  (Note  that  x increases  from  right  to  left  looking  from  below 
the  x-axls.)  The  cuts  along  which  the  graph  is  taken  are  shown  on  the 
rectangle.  In  the  case  of  multiple  cuts  in  the  same  direction,  a cut 
and  the  graph  along  the  cut  are  coordinated  by  the  same  line  type,  i.e., 
solid,  dashed,  etc. 
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4.3.2  The  ++  inodes 

Figures  4.2  through  4.6  are  plots  of  selected  ++  symmetry  natural 
modes.  Some  of  the  features  of  the  currents  in  Figure  4.2  are  common 
to  currents  of  all  symmetries.  The  mode  represented,  the  -H-l 
v/L  - 1.0  mode,  is  a first  resonance  as  manifested  by  the  single  "half- 
cycle" behavior  in  the  x-coordinate  cut  of  the  component  of  current. 
The  transverse  cut  - that  in  the  y-coordinate  - Indicates  a steeply 
peaked  behavior  of  the  current  near  the  edge.  This  is  the  representation 
obtained  by  the  graphical  interpolation  through  the  subsect ionally 
constant  method  of  moments  current  solution.  It  is  an  approximate 
numerical  representation  of  the  correct  singular  current  behavior  as 
given  by  (3.3).  The  method  of  moments  solution  also  manifests  an 
"undershoot"  in  the  first  current  zone  Inward  from  the  edge  zone  as 
indicated  near  y/L  « ±0.3  in  the  figure.  Because  of  the  rapid  variation 
of  this  hump  with  respect  to  wavelength,  we  hold  it  suspect.  A likely 
cause  for  such  anomalous  behavior  is  the  approximation  of  a singular 
function  by  a constant  zone  for  the  current  at  the  edge.  In  spite  of 
these  two  anomalies,  we  expect  the  subsectionally  zoned  current  to 
provide  good  approximate  scattered  field  Information  at  distances  from 
the  plate  on  the  order  of  nominally  one-half  wavelength  and  greater. 

Figures  4.3  through  4.5  give  the  modes  for  three  different  aspect 
ratios  on  the  ++2  pole  trajectory.  A comparison  of  Figure  4.4  with 
4.3  and  4.5  indicates  that  the  transverse  - J - current  is  much  larger 
in  magnitude  for  w/L  - 0.6  than  for  1.0  or  0.2.  Indeed,  it  is  the 
largest  transverse  component  along  the  -H-2  trajectory.  This  observation 
along  with  the  observed  loop  in  the  ++2  trajectory  is  Interpreted  to  be 


the  natural  mode  distribution 


natural  mode  distribution  for  the  ++2  mode  with  w/L 
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something  of  a "dual-resonance."  Namely,  the  aspect  ratio  w/L  - 0.6 
admits  to  a simultaneous  resonance  in  the  transverse  and  longitudinal 
current  components. 

The  common  feature  of  a three  "half-cycle"  longitudinal  current 
distribution  is  observed  in  all  three  -H-2  graphs.  This  is  consistent 
with  the  concept  of  a second  symmetric  resonance.  We  observe  in 
Figure  4.5  that  for  w/L  ■ 0.2  a breakup  in  the  transverse  current  solution 
is  manifested  by  rapidly  varying  phase  and,  generally,  oscillatory  be- 
havior of  Jy  w.r.t.  x.  This  is  interpreted  as  the  onset  of  numerical 
noise  for  the  narrow  strip  where  the  matrix  equations  become  ill- 
conditioned.  Recall  that  we  indicate  in  Section  3.4.2  that  the  maxi- 
mum transverse  zoning  density  which  produces  stable  and  consistent  pole 
results  is  used  -4x3  zones  in  this  case.  Because  the  transverse 
current  component  is  small  compared  with  the  longitudinal  current  com- 
ponent and  because  the  oscillations  in  the  computed  values  approximately 
cancel  under  integration,  we  expect  the  scattered  fields,  which  depend 
almost  solely  on  for  thin  strips,  to  be  acceptable  when  computed 
from  this  distribution. 

Figure  4.6  shows  the  ++3  w/L  - 1.0  natural  mode.  It  is  Included 
because  it  is  representative  of  modes  associated  with  the  -H-3  and  ++4 
trajectories.  The  characteristics  of  this  mode  are  similar  to  the 
++1  mode  except  that  the  transverse  component  of  current  is  larger  than 
the  longitudinal  component. 

4.3.3  The  +-  modes 

Figures  4.7  and  4.8  give  modes  associated  with  the  +-1  trajectory 
for  w/L  ■ 1.0  and  w/L  ■ 0.2,  respectively.  They  support  the  observation 
stated  in  Section  4.2  that  for  the  dominant  resonances  the  +-  modes 
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constitute  the  intermediate  resonances  between  the  ++  modes.  Comparisons 
of  Figure  4.7  with  4.2  and  4.3  and  of  Figure  4.8  with  4.5  show  a pro- 
gression from  one  "half-cycle"  on  ++1  to  two  on  +-1  and  then  three 
"half-cycles"  on  -H-2.  An  additional  feature,  which  is  common  to  +- 
and  — I-  modes,  is  observable  from  Figure  4.7.  Namely,  on  a square  struc- 
ture a +-  mode  or  a -+  mode  replicates  itself  under  a 90  degree  ro- 
tation. 


4.3.4  The  -+  modes 

As  Section  4.2  indicates,  the  —4-  and  — modes  are  related  through 
their  antisymmetries  in  the  narrow  direction  of  the  plate.  This  symmetry 
constraint  produces  a numerically  unstable  condition  in  the  narrow 
strip  limit. 

Figures  4.9  and  4.10  give  modal  distributions  for  the  — (-1  trajectory 
and  w/L  = 1.0  and  0.3.  The  features  of  modes  on  the  -+2  and  -+3  tra- 
jectories are  essentially  the  same.  The  — +1  w/L  = 1.0  curves  in  Figure  4.9 
indicate  antisymmetry  in  the  direction  transverse  to  current  flow.  In 
other  respects  the  mode  is  not  unlike  the  modes  discussed  in  previous 
subsections.  The  w/L  = 0.3  mode  indicates  a new  phenomenon,  however. 

The  breakup  in  the  mode  solution  for  narrower  strips  is  manifested  in 
the  dominant  current  component,  indicated  here  by  the  oscillations 
in  magnitude  and  phase  in  the  transverse  direction. 

Figure  4.11  gives  the  natural  mode  for  the  -+4  w/L  = 1.0  case. 

This  mode  shows  a three  "half-cycle"  variation  of  each  current  component 


in  the  direction  of  current  flow. 


asoMj 


4.3.5  The  — modes 


That  the  — n modes  are  Identical  to  the  ++n  modes  under  90  degree 

rotation  for  w/L  » 1.0  means  that  the  data  in  Figures  4.2,  4.3  and 

4.6  may  be  interpreted  as  — mode  data  with  the  roles  of  J and  J 

x y 

interchanged.  As  a companion  to  Figure  4.2,  we  present  the  — 1 w/L  = 0.2 
mode  in  Figure  4.12.  It  is  indicated  in  several  previous  sections  that 
for  narrower  strips  the  — solution  becomes  unstable.  This  instability 
is  evident  in  the  unusual  curve  shapes  indicating  a depletion  of  cur- 
rent near  the  center  of  the  strip  and  in  the  rapid  variation  of 
w.r.t.  x. 


Figures  4.13  and  4.14  show  modes  on  the  — 6 trajectory.  It  is 
noteworthy  that  for  w/L  **  1,  the  component  of  current  dominates. 

The  component  is  small  and  manifests  a rapid  variation  in  phase 
w.r.t.  y.  The  w/L  « 0.3  mode  on  this  trajectory  indicates  the  expected 
breakup  in  a subtle  fashion  only.  That  is,  the  component  of  current 
does  not  tend  toward  a singularity  at  x/L  « ±0.5  as  expected  and  as 
other  cases  manifest. 


4.4  Conclusions 

The  poles  and  specimen  natural  modes  for  the  dominant  resonant 
region  of  the  rectangular  aperture  are  given  in  this  chapter  . The 
computed  pole  trajectories  lend  themselves  to  physically  consistent 
interpretation.  The  broadband  character  of  square  and  near-square 
structures  noted  by  other  workers  is  seen  to  be  due  to  the  relative 
proximity  of  several  poles  to  the  imaginary  axis  in  the  complex  fre- 
quency plane.  Pole  trajectories  for  modes  that  have  a dominant  current 
component  that  is  antisymmetric  transverse  to  current  flow  are 


natural  mode  distribution  for  the  — 1 mode  with 


natural  mode  distribution  for  the  — 6 mode  with  w/L 
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undeterminable  by  the  present  method  for  strips  narrower  than  about 
w/L  - 0.3.  Because  of  their  antisymmetry,  they  are  expected  to  con- 
tribute negligibly  to  currents  In  problems  of  practical  Interest. 

The  natural  modes  are  characterized  here  by  specimen  data.  The 
modes  are  observed  generally  to  be  consistent  with  physical  expectations. 
Exceptions  to  this  statement  occur  for  the  unstable  case  referred  to 
above.  The  satisfactory  results  are  gratifying  In  light  of  the  low 
zoning  density  In  all  of  the  solutions.  The  data,  obtained  for  4x4 
zones /quadrant  by  and  large,  should  provide  good  scattered  field  esti- 
mates outside  a region  near  the  surface  of  the  structure.  The  resolution 
of  the  modes  can  be  improved  upon  by  using  the  pole  data  herein  as 
iteration  starting  points  from  which  to  conduct  higher-order  solutions. 

The  judgments  on  the  reliability  of  the  computed  mode  data  stated 
above  are  based  on  the  author's  experience  with  previous  method-of- 
moments  solutions.  The  conclusions  need  to  be  tested  for  a particular 
application  before  they  are  applied.  In  particular,  two  areas  must  be 
tested: 

1.  The  distance  by  which  an  observation  point  must  be  separated 
from  the  structure  In  order  to  give  smooth  scattered  fields 
from  the  coarsely  grained  current  expansion  must  be  determined; 
and 

2.  The  negligibility  of  — H and  — symmetry  modes  for  thin  structures 
must  be  tested.  This  may  be  done  by  observing  the  trend  in 
coupling  coefficients  for  these  modes  as  w/L  decreases. 
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5.  CONCLUSIONS 

We  attempt  to  fulfill  two  purposes  with  the  work  reported  here. 

The  first  Is  to  characterize  the  dual  problems  of  the  rectangular 
aperture  and  the  rectangular  plate.  The  second  Is  to  clarify  some  of 
the  previously  unsettled  issues  in  the  theory  of  SEM  regarding  suit- 
ability of  particular  coupling  coefficient  forms  and  the  role  of  the 
entire  function  portion  of  the  Mittag-Leffler  form  of  the  frequency 
domain  singularity  expansion. 

That  the  first  of  these  goals  is  met  is  supported  by  the  inclusion 
of  physically  plausible  SEM  data  given  herein.  Both  the  poles  and 
natural  modes  determined  lend  themselves  to  straightforward  physical 
interpretation.  The  natural  mode  information  suffers  from  a "graininess" 
due  to  the  use  of  minimal  zoning  density  in  the  numerical  procedures. 

It  is  anticipated  that  the  only  limitation  imposed  thereby  is  that 
coupled  or  scattered  fields  cannot  be  predicted  from  the  modes  in  the 
immediate  vicinity  of  the  body. 

The  scope  of  the  present  effort  did  not  allow  time  domain  expansion 
of  fields  from  the  data  herein.  The  ultimate  judgment  on  the  correct- 
ness and  applicability  of  these  results  rests  in  their  time  domain  ap- 
plications. This  testing  shall  necessarily  include  the  causality  tests 
proposed  in  Chapter  2. 

The  second  goal  has  been  partially  fulfilled.  Again,  ultimate  judg- 
ment must  be  deferred  until  the  theoretical  concepts  in  Chapter  2 are 
tested  for  their  completeness  in  a broader  range  of  problems.  Specifically, 
the  conclusions  herein  on  these  issues  may  be  summarized  as  follows. 
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1.  In  a numerically  conducted  SEM  solution  neither  the  entire 
function  nor  the  asymptotic  behavior  of  the  frequency  domain 
inverse  kernel  Is  explicitly  recoverable. 

2.  If  the  Laplace  transform  Inversion  procedure  Is  applied  to 
the  inverse  kernel  as  a unit,  assuming  the  net  asymptotic 
behavior  is  available,  the  entire  function  contributes,  at 
most,  impulsive  terms  at  the  "turn-on"  of  the  current. 

3.  The  effect  of  the  asymptotic  behavior  of  the  frequency  inverse 
kernel  is  to  dictate  the  support  of  its  time  domain  counterpart. 
Functionally,  in  SEM  this  affects  the  form  of  the  coupling  co- 
efficient. 

A.  The  class  1 coupling  coefficient  and  the  asymptotic  form  im- 
plicitly assumed  in  using  it  are  potentially  inconsistent  with 
causality  if  taken  alone. 

5.  We  have  developed  a conservative  "class  X"  coupling  coefficient 
which  imposes  causality  directly.  Further,  we  describe  numer- 
ical tests  whereby  one  may  ascertain  the  consistency  with 
causality  of  a given,  less-conservative  coupling  coefficient, 
such  as  the  class  1 or  class  2 forms. 

6.  The  asymptotic  behavior  of  the  series  portion  of  the  Mittag- 
Leffler  form  of  the  frequency  domain  expansion  may  be  modified 
to  admit  to  a class  1 coupling  coefficient  representation. 

This  can  quite  possibly  lead  to  acausal  contributions  to  the 
solution,  however,  and  an  entire  function  contribution  with 
support  over  a finite  time  must  be  included  in  the  total  solu- 
tion in  order  to  provide  results  consistent  with  causality. 
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We  note  several  times  in  Chapter  2 that  the  class  X coupling  coef- 
ficient may  prove  unduly  conservative  for  some  problems.  In  fact,  it 
is  conservative  for  the  transmission  line  and  linear  scatterer  examples 
presented.  Further  work  is  needed  to  determine  the  conditions  under 
which  the  class  X imposition  of  causality  may  be  relaxed  to  the  class  2 
form.  The  asymptotic  behavior  associated  with  the  class  2 coefficient  is 
equivalent  to  the  asymptotic  behavior  of  the  residue  series  term  by 
term.  However,  this  series  is  not  uniformly  convergent,  in  general, 
and  therefore  does  not  admit  directly  to  termwise  integration.  However, 
in  the  time  domain  inverse  kernel  plots  in  Chapter  2,  the  clear  time 
between  time  zero  and  the  time  at  which  causality  predicts  the  arrival 
of  an  incident  wave  indicates  that  the  termwise  integration  would, 
indeed,  produce  correct  results.  The  question  of  whether  or  not  there 
is  ever  a need  to  include  an  entire  function  contribution  with  the 
class  2 and  class  X coupling  coefficient  expansions  remains  an  open 


question. 
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